復小波域混合概率圖模型的超聲醫學圖像分割

夏平 施宇 雷幫軍 龔國強 胡蓉 師冬霞

醫學影像已成為醫生臨床疾病輔助診斷最重要的手段之一[1?3];超聲醫學影像以其實時顯示、對人體無傷害、適合軟組織診斷等特點在眼部、心臟、腹部以及人體淺表部位等部位的疾病診斷中廣泛應用.超聲醫學圖像分割是醫學影像分析的首要前提,是獲取醫學影像中病變區域與特定組織、器官信息的一種重要手段,其目的是從超聲影像中分離出醫生診斷感興趣的區域,對輔助醫生臨床診斷、定位病變組織、規劃手術治療等方面顯示出越來越重要的臨床價值[4].

由于超聲成像過程中的干涉現象及采集設備電子器件的隨機擾動,使得超聲圖像存在大量不規則的斑點噪聲,降低了圖像中目標與背景區域的對比度,大大降低了超聲圖像的分辨率,使得圖像中不同組織間的灰度差異性不明顯,削弱了超聲圖像質量,從而影響了醫生對病灶的識別,也增大了醫生查出彌散性組織病變的難度.獲取超聲醫學圖像中特定組織器官或病變區域定量信息的過程中,超聲醫學圖像分割是不可或缺的手段[5];超聲影像的分割除傳統分割手段[6]外,近年來隨著統計學理論[7]、模糊集理論[8]、神經網絡[9]、多尺度理論[10]的發展與應用出現了很多新的分割算法[11?14].文獻[15]針對醫學圖像分割中區域主動輪廓模型易受目標和背景區域面積比的影響,且對初始位置敏感的問題,提出模糊C-均值聚類(Fuzzy C-means,FCM) 協作改進CV (Chan-Vese) 模型(FCM-CV) 的圖像分割算法;該算法通過在CV 模型中增加能量權值函數消除面積比的影響后,用FCM 粗分割結果指導設定改進CV 模型零水平集的初始位置,從而提高了圖像分割的穩定性和準確性.隨著小波技術的發展與應用,基于小波域馬爾科夫隨機場(Markov random field,MRF) 分割方法的研究十分活躍;文獻[16]提出了小波域中樹結構化MRF 的圖像分割算法,按照圖像分類層次樹的結構形式,該算法在小波域多分辨率中構建了序列MRF 模型,它們相互嵌套,并通過在小波域的相鄰尺度間和同一分辨率內的兩個層次上遞歸運算實現圖像的有效分割;該算法在表示圖像結構信息的同時,較好地描述了圖像的非穩態性,得到較好的分割效果.但小波分析本身缺乏平移不變性,且只能在水平、垂直、以及45?三個方向建模圖像的特征信息,因而,基于小波的圖像分割在獲取目標邊緣的連續性不理想,檢測結果存在一定“空洞”現象.相對于小波分析,雙樹復小波變換[17]采用兩路實離散小波分析實現復小波分析,改善了小波分析中對平移的敏感性,并實現了信息的多方向獲取,在一定程度上克服了目標分割的“空洞”問題;文獻[18]利用雙樹復小波矩形狀信息對醫學淋巴結進行精確定位,對已探測的淋巴結進行形狀特征提取,判斷其是否為真實的淋巴結;然后,應用水平集DRLSE(Distance regularized level set evolution) 模型對淋巴結進行輪廓描述,實現淋巴結完整的分割.

超聲醫學影像中,其成像過程所造成的斑點噪聲以及相關組織紋理是影響圖像分割準確性的重要因素,結合超聲影像的特點,本文提出了復小波域混合概率圖模型的超聲醫學圖像分割方法.利用雙樹復小波變換(Dual tree-complex wavelet transform,DT-CWT) 的多分辨率分析具有的多方向性和近似的平移不變性,捕捉超聲圖像的局部弱特征信息[19];多分辨率的相鄰層間構建的有向圖模型確立父-子節點標號相互聯系,構建MRF 模型來建立同層復小波系數標號與鄰域間的聯系;而噪聲是統計獨立的,從而,通過構建混合概率圖模型建立了當前節點的標號與其祖先節點及同層鄰域間節點之間聯系,有效地獲取超聲影像中有用信息,實現超聲影像的可靠分割.

1 超聲影像的雙樹復小波域稀疏描述

DT-CWT 采用兩路實小波對信號進行分解以實現復小波變換[17],其中,一路為CWT 的實部,另一路為其虛部;復小波函數:Ψ(t)=Ψh(t)+jΨg(t),其中,Ψh(t)、Ψg(t) 均為實小波.其分解過程用樹A和樹B 兩個實小波分解樹表述,對應樹A、樹B的濾波器分別為h0(k)、h1(k) 和g0(k)、g1(k);如圖1(a).每層CWT 分解得到2 個低頻分量和±15?、±45?、±75?等6 個高頻方向子圖像.如圖1(b) 所示.

樹A、樹B 對應的濾波器組除滿足完全重構的條件外,還需使Ψ(t) 解析[17],即:

將樹A、樹B 兩個實DWT 分別由兩個方陣Hh和Hg表示,則DT-CWT:

超聲圖像I進行復小波分析,有:Ch=Hh·I和Cg=Hg·I;其中,Ch、Cg分別為DT-CWT 的實部和虛部系數,Ch+jCg為DT-CWT 稀疏描述的高頻復值方向子帶,從而得到6 個方向子帶:

式中,每一子帶的復小波實部與虛部由v={1,2}確定,mmm=(m1,m2) 為位置坐標,6 個方向子帶用u ∈{1,2},i ∈{1,2,3} 標示.

圖1 DT-CWT 變換及其子帶方向Fig.1 DT-CWT and the sub-band direction

相對于實小波分析,DT-CWT 繼承了小波分析的多尺度分析特性和稀疏表示能力,同時具有更好的方向選擇性與近似的平移不變性.

2 復小波域超聲圖像的混合概率圖模型

N × N格網S上的超聲醫學圖像I,經J -1 級復小波分析后形成J幅不同尺度的高頻子帶圖像,J幅圖像大小集合為其中,S(0)為原始圖像;S(j)對應第j級復小波分析的子帶圖像.本文基于格網集合上構建復小波域超聲醫學圖像的混合概率圖模型.如圖2 所示.

2.1 標記場的混合概率圖模型建模

定義復小波系數的標記場:

式中,標記FS(n)對應于分辨率n的復小波系數 CS(n) 的標記場,標記場F(n)=表示尺度n上空間(i,j) 的標記值.類似地,FFFS(J?1)對應最低分辨率復小波系數CCCS(J?1)的標記場,FFFS(0)對應于原始圖像的標記場.

標記場描述了像素點特征值在“父-子”間的屬性傳遞性以及同分辨率中局部的相關屬性,不同分辨率尺度上的觀測特征展現對應的區域特征,因而標記場的多分辨率描述體現了這兩種屬性共同作用的結果.

依據復小波域系數的持續性,相鄰尺度間復小波系數的標記具有Markov 性,即多分辨率分析構成的層間Markov 鏈表征了層間的相關性,并且,將父節點標記屬性傳遞給其子節點,尺度間的“父-子”節點聯系形成復小波域中系數間的貝葉斯網絡結構.如圖3(a) 和(b) 所示.

圖2 DT-CWT 域系數向量結構Fig.2 Coefficient vector structure of DT-CWT domain

標記場建模中,尺度間標記具有一階Markov轉移特性[20].多分辨率分析中,相鄰尺度間的標記場服從一階Markov 分布,用轉移概率描述父、子復小波系數標號的依賴關系為:

式(5) 的一階Markov 性反映了高分辨率尺度上的標記場僅受相鄰的低分辨率尺度標記場影響,與其他尺度的標號無關.如圖3(c) 所示.

圖3 DT-CWT 域標記場貝葉斯網絡模型Fig.3 Marking field Bayesian network model in DT-CWT domain

尺度間標記的傳遞屬性僅考慮相鄰層間“父-子”節點的相互作用,如圖3(b) 中父節點“a”與子節點“b”等;節點標記傳遞屬性數學建模:

式中,fi+1、fi分別表示父、子節點標記值.設醫學圖像分為K類,其傳遞屬性用父節點標記值到其子節點標記值轉移概率p來描述,則標記轉移概率矩陣:

式中,ε ∈[0,1]為相鄰尺度間“父-子”節點具有相同標記的概率.

式(8) 轉移概率定義使式(7) 構成一強對稱矩陣,此時,“父-子”節點標記在傳遞過程中信息容量達到最大.本文采用ε=1,即高分辨率尺度標記場的初始值由相鄰低分辨率尺度標記場的實現來確定.

標記場建模中,尺度內標記具有Markov 性.即尺度n中復小波系數的標記場是以?ij為鄰域的MRF,即:

式中,Pf(n)>0,?ij為尺度n內的(i,j) 位置的鄰域系統.

同一尺度標記場是以鄰域系統的結構作為條件,MRF 建模中常采用1 階或2 階鄰域系統,本文采用2 階鄰域系統描述同一尺度標記場的結構間聯系.如圖4 所示.文獻[20]已證明,鄰域系統基團勢函數參數αm、βi、γi、δ等與標記場取值無關,因此,本文在復小波域MRF 建模中,對指定階數的所有鄰域系統均采用相同的勢函數參數.

圖4 DT-CWT 域標記場的2 階鄰域系統Fig.4 The second order neighborhood system in DT-CW

據Hammersley-Clifford 定理[21],MRF 和Gibbs 分布具有等效性,使用MLL (Multilevel logistic model) 模型來表示標記場的先驗分布[7,12].本文采用各向同性的雙點MLL 模型描述像素空間鄰域的作用關系,僅考慮其二元勢函數來計算標記之間的局部約束條件.

本文中,定義隨機場?={?ij|(i,j)∈S=N×N}是S的鄰域系統,?ij是(i,j) 的鄰域點集,用各向同性的雙點MLL 模型來計算標記之間的局部約束條件,其勢函數Uc與標記場的局部條件概率P分別用式(10) 和式(11) 描述.

式中,Z為概率分布歸一化因子,即分配函數,c是二維圖像基團,β是對應的勢團參數,δ(·,·) 為信息采樣函數,具有如下特性:是位置(i,j) 的鄰域中不等于的相鄰位置的數目.

尺度間與尺度內兩種標記場體現了多分辨率分析中不同層間、同一層內的復小波系數特征相互作用與影響后表現出的特征屬性,尺度間作用構成了一個有向圖結構,在此結構中,父節點標記屬性是以式(7) 轉移概率矩陣形式傳遞給子節點,因此,多層傳遞中形成了一個貝葉斯網絡結構:尺度內聯系構成了無向MRF 的圖模型,每個節點標記與鄰域?ij節點標記建立聯系,從而,在尺度內構成了無向圖模型:二者共同構建了復小波域空間標記場的混合概率圖模型.

2.2 觀測特征場建模

超聲圖像I做J-1 層復小波分解后高頻系數為CCC,每一尺度復小波用對應層號n(1≤n ≤J -1)表示.將每一分辨率上對應位置不同頻帶方向的復小波系數構成向量,形成了對應的向量圖像,則此向量圖像表示了該分辨率中某一方向圖像的觀測特征,層號n=0 對應原始圖像.則:

分辨率尺度n、位置(i,j) 的復小波系數的特征向量:

經J -1 層復小波分解后,定義圖像I的格網S結構:

由式(15) 可知,J -1 分解后格網中共包含6J+1 個子帶,尺度上的格網:

各尺度上對應的復小波系數C的向量場描述的特征場的位置在(i,j)∈S(n)所表示的范圍內.不同位置集上相應的復小波系數形成觀測特征場表示為:

式中,C(J?1)為最粗尺度的復小波系數,只有一個節點,C(0)對應原始圖像.

特征場建模中,考慮到任意節點上復小波系數值的概率由每一觀測特征由所屬標號和所在分辨率的觀測特征值確定,不同標號下觀測特征相互統計獨立.即:

圖5 給出了前列腺囊腫圖像復小波分析后第1、第2 層各方向高頻圖像(圖5(a)) 及系數統計圖(圖5(b)),統計分析表明,信號經DT-CWT 后子帶系數的邊緣概率分布呈現尖峰、長拖尾的非高斯分布,如圖5 所示,反映了DT-CWT 分析在尺度內對信號的“聚集”能力和稀疏表示能力.由此可知,子帶高頻系數值大部分位于“0”值附近,少數系數值較大,為此,用值“大”或“小”兩狀態描述子帶系數的分布規律;基于此,文獻[22]采用兩個不同方差的高斯混合模型來近似逼近子帶系數總體分布.本文中尺度n的觀測特征值采用高斯混合模型建模,該分辨率中每一標號的觀測特征值采用高斯模型擬合,則:

圖5 醫學圖像的DT-CWT 分解及其高頻系統統計Fig.5 Medical image of DT-CWT and Statistical characteristics of the coefficient

式中,分別表示尺度n上標號的均值向量和協方差矩陣;Q為特征向量維數.

2.3 誤分割率最小的超聲醫學圖像分割

分割的目標是對復小波域中每個復小波系數觀測特征與其分配的標記一一對應,記:(Cs,Fs);用條件概率P(Cs|Fs) 表示復小波系數對標記的依賴關系,用邊緣分布P(Fs) 表示標號的先驗概率分布,在給定復小波系數Cs下標記Fs的概率:

分割過程中標記的獲取與其觀測特征應相適應,其誤分割率Pe表示為:

最佳分割的目標是分配一組最佳標記Fs使圖像的誤分割率達到最小,即在已知觀測復小波系數Cs的基礎上,求使得P(Fs|Cs) 最大的分割

式(22) 等價于最小化能量問題.能量函數為標記場能量U和特征場能量U之和.即:

由式(20)～(24),采用基于局部條件概率的確定性算法,通過逐點更新圖像標記完成圖像分割.

2.4 分割算法實現

超聲圖像I大小為N ×N,對其進行J -1 層復小波分解,采用高斯混合模型建模各尺度特征場,每一標號的觀測特征場采用高斯模型擬合,用MLL模型建模各個尺度的標記場,使用EM 算法完成模型中參數估計.分割算法步驟:

1) 對超聲圖像進行復小波多分辨率分解,獲得觀測場,獲取最粗尺度上的初始分割.本文采用Kmean 聚類算法獲取尺度n=J -1 上的初始分割結果;

2) 令尺度n ←n-1;

3)E步.用MPL 方法估計式(25)、式(26) 的模型參數[8]:

式中,n表示多分辨率分解的尺度,l表示迭代次數,m表示標記類別;

4)M步.用式(25) 和式(26) 估計的參數,采用ICM 求解式(23) 的最小能量,獲取該尺度上標記場結果;

5) 尺度內重復步驟(3) 和步驟(4),直到滿足迭代條件(?U=Ui-Ui?1≤10?3) 即停止迭代,得尺度n的標記場;

6) 尺度間迭代:將尺度n的標記場結果映射到最鄰近的精細尺度n-1 上,作為尺度n-1 的初始分割,重復步驟2),直至獲取最細尺度的標記場;

7) 進行形態學處理,實現超聲醫學圖像分割.

3 實驗仿真與性能分析

3.1 實驗參數與定量評價指標

采用囊腫、乳腺癌、輸尿管末端囊腫、左腎腫物、前列腺等超聲醫學圖像進行實驗,醫學圖像源由市中心人民醫院影像科提供.對比實驗時采用復小波域K-means 算法、空域MRF 算法、小波域MRF 算法以及本文算法對圖像進行分割.

參數選取:復小波分析采用Q-shift DT-CWT,進行3 層DT-CWT 分解;K-means 算法中,分為目標與背景2 類;MRF 中,勢團參數選β=0.8;迭代運算條件:前后兩次迭代能量之差?U ≤10?3.

定量評價指標采用概率Rand 指數(Probabilistic rand index,PRI)[23]、信息變化指數(Variation of information,VoI)[24]、全局一致性誤差指數(Global consistency error,GCE)[25]、邊界偏移誤差指數(Boundary displacement error,BDE)[26]等客觀評價指標,從不同的角度反映分割結果的優劣以及對算法的比較與評價.

1) 概率Rand 指數(PRI) 用于度量待評測結果與參考結果之間的屬性共生的一致性.待評測聚類標記場FFF,圖像中的每個像素點i標記號fi,用信息采樣函數δ(·,·) 指示兩個像素點i、j是否屬同一類.

參考分割結果記為?,用集合?={r11,r12,···,rMN} 描述,則待評測分割結果與參考結果之間PRI定義[23]:

PRI ∈[0,1],PRI值越大兩個聚類結果的一致性越高,分割結果越好.

2)V oI從信息論的角度度量不同聚類之間的距離,反映分割信息量的丟失.定義[24]:

式中,I(F,?) 是F和?所共有的互信息.V oI ∈[0,∞),V oI值越小兩個聚類結果的一致性越高,分割效果越好.

3)GCE用于度量不同分割結果之間的一致性誤差[25].GCE ∈[0,1],GCE值越小表示細化誤差越小,分割效果越好.

4)BDE待評測結果和參考結果邊緣像素間平均距離來度量分割結果[26].BDE ∈[0,∞),值越小表示兩者邊界差異越小,分割效果越好.

3.2 實驗結果與性能分析

結果如圖6、圖7 所示.由圖6、圖7,從清晰度角度看,基于空域MRF、小波域MRF 以及本文算法分割的圖像均較清晰;從分割區域效果看,本文算法很好地將目標區域分割出來,其次,表現較好的是小波域MRF 算法,而其他算法表現差強人意;其主要原因在于,復小波域K-means 分割中并沒有很好地利用復小波分析優良的特征提取特性;而空域MRF 算法中,在圖像分割區域內部構建MRF,充分考慮了區域內部的一致性問題,但忽視了不同目標區域邊界的各向異性,從而,導致應用該算法進行圖像分割時易產生目標邊緣模糊;小波域MRF 算法相比于本文算法,主要缺乏平移不變性以及方向特征信息表示不足.4 種分割算法的定量比較見圖6,通過PRI、V oI、GCE、及BDE等指標的評價,本文算法的PRI指標高于其他算法0.014 以上;V oI指標低0.09 以上;GCE指標低0.03 以上;BDE指標低1.10 以上.本文算法的分割結果從邊緣、輪廓的清晰度、細節信息捕捉等方面均好于其他三種算法.

圖6 超聲醫學圖像4 種算法分割結果比較Fig.6 Comparison of segmentation results of four algorithms for ultrasonic medical images

圖7 醫學超聲圖像4 種算法定量評價指標比較Fig.7 Comparison of quantitative evaluation indexes of four algorithms for medical ultrasound images

由圖6 和圖7,從清晰度角度看,基于空域MRF、小波域MRF 以及本文算法分割的圖像均較清晰;從分割區域效果看,本文算法很好地將目標區域分割出來,其次,表現較好的是小波域MRF算法,而其他算法表現差強人意;其主要原因在于,復小波域K-means 分割中并沒有很好地利用復小波分析優良的特征提取特性;而空域MRF 算法中,在圖像分割區域內部構建MRF,充分考慮了區域內部的一致性問題,但忽視了不同目標區域邊界的各向異性,從而,導致應用該算法進行圖像分割時易產生目標邊緣模糊;小波域MRF 算法相比于本文算法,主要缺乏平移不變性以及方向特征信息表示不足.4 種分割算法的定量比較見圖6,通過PRI、V oI、GCE、及BDE等指標的評價,本文算法的PRI指標高于其他算法0.014 以上;V oI指標低0.09 以上;GCE指標低0.03 以上;BDE指標低1.10 以上.本文算法的分割結果從邊緣、輪廓的清晰度、細節信息捕捉等方面均好于其他三種算法.

4 結論

本文提出了復小波域多分辨率MRF 模型的超聲醫學圖像分割算法.復小波分析中每一分辨率提供了6 個方向的高頻特征信息,并具有近似的平移不變性,適應了圖像的非平穩性,從而能較好地捕捉超聲圖像弱特征信息;誤分割率最小算法充分考慮復小波域中特征信息的層間、層內信息的相關性,利用了復小波域中每一尺度內同標號區域的特征信息分布規律,彌補了空域MRF 分割中參數估計復雜、小波域MRF 分割中缺乏平移不變性和特征提取方向性差的不足,得到了更好的分割效果.

理論推導和實驗仿真驗證了本文方法應用于超聲醫學圖像目標檢測的有效性以及實現超聲醫學圖像目標檢測所具有良好的檢測性能.

未來及目前正著手做的工作包括:1) 本文算法在實現誤分割率最小的能量函數最優解的求解中采用ICM 算法,正著手引入圖割、置信度傳播等算法實現能量最優解,探討通過優化求解能量最小化模型,提高算法的分割效率;2) 探討其他算法,如深度學習方法,應用于細胞、視網膜等醫學圖像分割;3)針對肺部CT 影像,在三維測量基礎上,結合本文算法研究檢測肺結節的方法預測和分析肺部腫瘤惡性情況,為臨床診斷與治療提供輔助手段.