本發(fā)明涉及頸動脈超聲圖像處理領(lǐng)域，尤其是頸動脈超聲圖像的內(nèi)中膜厚度測量算法。具體講,涉及頸動脈超聲圖像內(nèi)中膜厚度測量方法及系統(tǒng)。

背景技術(shù)：

超聲成像中頸動脈血管壁的內(nèi)中膜厚度(Intima Media Thickness,IMT)可以作為評估心血管疾病的早期病變程度的重要指標(biāo)，對突發(fā)心肌梗塞、中風(fēng)的診斷和預(yù)防有很重要的價值。傳統(tǒng)上，頸動脈超聲圖像的內(nèi)中膜厚度通常由人工手動測量獲得，測量者手工描繪圖像中的管腔-內(nèi)膜邊界(Lumen-Intima Interface,LII)和中膜-外膜邊界(Media-Adventitia Interface,MAI)，然后通過計算得出兩邊界間的距離獲得IMT。然而這種測量方式工作量大，非常耗時并且最終獲得的結(jié)果受測試者所受訓(xùn)練、個人經(jīng)歷及其主觀判斷的影響。手動測量的缺陷促進(jìn)了計算機輔助的IMT測量方法的發(fā)展。本發(fā)明提出了一種結(jié)合水平集分割和馬爾可夫隨機場(Markov Random Field,MRF)模型的頸動脈超聲圖像內(nèi)中膜厚度測量算法。

基于幾何變形模型的水平集方法由Osher和Sethian于1982年提出。該方法是利用偏微分方程作為數(shù)值分析放過與技術(shù)手段被廣泛運用于輪廓面或輪廓線的運動跟蹤。該方法將低維閉合曲線演化的問題轉(zhuǎn)化為高維空間中誰破解函數(shù)全面演化的隱含方式來求解，因此可以適應(yīng)于拓?fù)浣Y(jié)構(gòu)變化的處理。水平集方法的思想中最重要的就是Hamilton-Jacobi方程，它為運動的隱式曲面進(jìn)行基于時間的方程的數(shù)值求解。該方法的基本原理是，將曲線或曲面作為零水平集嵌入到高一維的水平集函數(shù)中，通過一個更高維的函數(shù)來表達(dá)低維的曲線或曲面的演化過程。零水平集被嵌入到水平集函數(shù)中，只要確定零水平集位置，就可以確定運動曲線或曲面演化的結(jié)果。下面簡要分析曲線模型的演化過程。平面閉合曲線C(s,t)＝(m(s,t),n(s,t)):0≤s≤1(s為弧長參數(shù)，t表示時間，m、n是平面曲線的坐標(biāo))被水平集方法隱含地表達(dá)為三維連續(xù)函數(shù)曲面φ(m,n,t):的一個具有相同函數(shù)值的等值曲線，令函數(shù)φ為符號距離函數(shù)，可將曲線C的變化歸結(jié)為由于函數(shù)發(fā)生了某種相應(yīng)的變化所致。因此，可將隨時間變化的閉合曲線表示為水平集函數(shù)對時間的變化：

上式把對曲線的演化轉(zhuǎn)換為對水平集函數(shù)的分析。

可用如下的偏微分方程來表示曲線C沿法線方向的演化過程：

上式被稱為水平集方程(Level Set Equation)，函數(shù)V被稱為速度函數(shù)(Speed Function)，表示曲線上個點的演化速度，N是曲線的法線方向。對式(2)求取全微分，得：

式中是φ的梯度，與曲線的法向方向相同。假設(shè)Ω為整幅圖像，水平集函數(shù)位于曲線C以內(nèi)的部分為負(fù)，以外的部分為正，則水平集的內(nèi)向單位法向量為

將式(2)和式(4)代入式(3)中，整理可得

這樣就用水平集方程表示了曲線演化過程，即曲線演化方程的歐氏表達(dá)，它是一種Hamilton-Jacobi類型的偏微分方程。給定幾何活動輪廓模型的偏微分方程，以及初始的水平集函數(shù)φ(m,n):φ(C)＝0，式(5)就可以保證水平集函數(shù)φ(m,n,t)隨時間的演化滿足{φ(C(s,t),t)＝0}條件，即φ(m,n,t)的零水平集始終是輪廓線C(s,t)。

技術(shù)實現(xiàn)要素：

為克服現(xiàn)有技術(shù)的不足，本發(fā)明旨在提出一種可靠、自動頸動脈超聲圖像的測量方法，使用計算機輔助測量的方式避免手動測量方式存在的缺陷。本發(fā)明采用的技術(shù)方案是，頸動脈超聲圖像內(nèi)中膜厚度測量方法，步驟如下：

1)圖像裁剪，裁剪移除無關(guān)信息，只保留血管部分；

2)自動提取感興趣區(qū)域ROI(Region of Interest)；

3)圖像濾波，濾除回波反射帶來的斑塊噪聲；

4)水平集分割，采用水平集分割算法追蹤血管壁的管腔-內(nèi)膜邊界LII(Lumen-Intima Interface,)；

5)初始標(biāo)記場，根據(jù)頸動脈血管壁各層近似平行的特性，在ROI中，將血管壁內(nèi)膜邊界向下平移一定像素距離，作為中膜-外膜邊界中膜-外膜邊界MAI(Media-Adventitia Interface,)的初始狀態(tài)，然后依據(jù)這兩條邊界將ROI分割為管腔、內(nèi)中膜、外膜三個部分，并將獲得的分割圖像作為Markov模型的初始標(biāo)號場；

6)再分割，引入Markov模型，根據(jù)初始標(biāo)號場和灰度場對圖像再次分割圖像；

7)后處理，對分割圖像后處理，移除錯分小區(qū)域，獲得最終分割圖像和LII、MAI，并測量內(nèi)中膜厚度IMT(Intima Media Thickness)。

步驟2)具體步驟是：根據(jù)超聲圖像灰度的統(tǒng)計規(guī)律，采用最大類間差法自適應(yīng)獲取灰度閾值，二值化超聲圖像；然后，對二值化超聲圖像形態(tài)學(xué)閉運算，再填補圖像中的孔洞，利用連通域面積較大以及質(zhì)心縱坐標(biāo)較大的原則，選擇期望的白色區(qū)域，其與黑色區(qū)域相交的上邊界認(rèn)為是LII的近似位置；最后以LII的近似位置為標(biāo)準(zhǔn)，向上取40、向下取20像素尺度，分別作為ROI的上、下邊界，獲得兩個分別對應(yīng)原始圖像ROI。

步驟3)圖像濾波是基于像素間距離相似性和灰度相似性的濾波，濾波輸出為

其中g(shù)(r)是濾波后圖像，是像素ξ和像素r的相似性函數(shù)，為歸一化函數(shù)，σ₁是定義域方差，σ₂是值域方差；

步驟4)進(jìn)一步及細(xì)化為：

Step1：初始化，假設(shè)t_i時刻的水平集函數(shù)φ(m,n,t_i)為已知，m、n是平面曲線的坐標(biāo)；

Step2：求解水平集方程：獲得t_i+1時刻的水平集函數(shù)φ(m,n,t_i+1)在整個圖像區(qū)域的值，此時記為Γ_i+1的演化曲線就是φ(m,n,t_i+1)的零等值曲線，即Γ_i+1＝{φ(m,n,t_i+1)＝0}，此時的φ(m,n,t_i+1)不再是符號距離函數(shù)；

Step3：重新初始化，使用φ(m,n,t_i+1)代替下式中的初始水平集函數(shù)φ₀，迭代求解至穩(wěn)定解，仍記為φ(m,n,t_i+1)：

其中φ_t為t時刻的水平集函數(shù)，φ(m,n)為初始的水平集函數(shù)；

Step4：結(jié)合φ(m,n,t_i+1)的值，求解物理量的控制方程，得到t_i+1時刻的物理量的值；

Step5：重復(fù)步驟，進(jìn)入下一時刻的計算，直至停止。

步驟6)建立待分割圖像的Markov模型，對于Markov模型的求解，流程如下：

Step1:輸入待分割的原始圖像；

Step2:對圖像進(jìn)行初始化，得到圖像的初始標(biāo)記場；

Step3:在上一步標(biāo)號場的基礎(chǔ)上對圖像中每個像素遍歷求使之取到滿足迭代條件的標(biāo)記，更新原有的標(biāo)號場；

對Step3進(jìn)行迭代直到收斂條件。

步驟6)進(jìn)一步細(xì)化為：

圖像分割是基于像素的特征屬性和區(qū)域?qū)傩越o每一個像素分配標(biāo)號的過程。在馬爾可夫隨機場模型(Markov Random Field,MRF)中，用兩個隨機場來描述待分割的圖像，一個是標(biāo)號場X，常稱為隱隨機場，用先驗分布P(X)描述標(biāo)號場的局部相關(guān)性；另一個是灰度場或特征場Y，以標(biāo)號場為條件，用分布函數(shù)P(Y|X)描述觀測數(shù)據(jù)或特征向量的分布，其中標(biāo)號場X是馬爾可夫場，滿足Markov特性，其先驗分布表示為

其中A＝{1,…a}是圖像中像素的集合，a是圖中像素總數(shù)，x_i∈{1,…,q}是像素i的類別標(biāo)號，q是圖像中像素類別的總數(shù)，N_i為像素i的鄰域，根據(jù)Hammersley-Clifford等價定理，所有像素的聯(lián)合概率分布為

其中Z是歸一化常數(shù)，U(x)是先驗?zāi)芰?/p>

c是一個基團，即包含若干位置的集合，x_i、x_j是像素點i、j的類別標(biāo)號，極端情況下認(rèn)為鄰域結(jié)構(gòu)的所有子集都是基團，C表示基團的集合，V_c(x)是定義在基團c上的勢函數(shù)，它僅與像素的鄰域像素有關(guān)，式中V₁(x_i)對應(yīng)于鄰域相關(guān)函數(shù)，勢團模型選擇二階的Potts模型進(jìn)行估計，即僅考慮二階鄰域系統(tǒng)中兩點間的相互作用，則用偽似然近似逼近的先驗概率表達(dá)為：

Potts模型：

其中β是控制鄰域間作用強度的函數(shù)，N_i、n_i是像素i的鄰域像素；

另一個隨機場Y指待分割圖像，它是可觀測的隨機過程，即P(Y)是固定的，則根據(jù)貝葉斯定理：

P(X|Y)∝P(X)P(Y|X) (13)

將模型轉(zhuǎn)化為對似然函數(shù)P(Y|X)和先驗概率P(X)的求解，這里P(X|Y)是圖像的后驗概率，其中條件分布函數(shù)P(Y|X)，通常用高斯分布描述：

其中為相同標(biāo)記為x_i的區(qū)域像素灰度的方差，為像素i鄰域灰度的均值；

引入Markov模型，根據(jù)公式(11)、公式(14)求解模型。

選取Potts模型作為公式(10)中的勢團表達(dá)函數(shù)，采用條件迭代算法ICM(Iterative Conditional Mode)求解Markov模型。

頸動脈超聲圖像內(nèi)中膜厚度測量系統(tǒng)，由超聲波裝置及計算機組成，超聲波裝置測得的圖像信息發(fā)送給計算機處理，計算機設(shè)置有如下模塊：

1)圖像裁剪模塊，裁剪移除無關(guān)信息，只保留血管部分；

2)自動提取感興趣區(qū)域ROI(Region of Interest)模塊；

3)圖像濾波模塊，濾除回波反射帶來的斑塊噪聲；

4)水平集分割模塊，采用水平集分割算法追蹤血管壁的管腔-內(nèi)膜邊界LII(Lumen-Intima Interface,)；

5)初始標(biāo)記場模塊，根據(jù)頸動脈血管壁各層近似平行的特性，在ROI中，將血管壁內(nèi)膜邊界向下平移一定像素距離，作為中膜-外膜邊界中膜-外膜邊界MAI(Media-Adventitia Interface,)的初始狀態(tài)，然后依據(jù)這兩條邊界將ROI分割為管腔、內(nèi)中膜、外膜三個部分，并將獲得的分割圖像作為Markov模型的初始標(biāo)號場；

6)再分割模塊，引入Markov模型，根據(jù)初始標(biāo)號場和灰度場對圖像再次分割圖像；

7)后處理模塊，對分割圖像后處理，移除錯分小區(qū)域，獲得最終分割圖像和LII、MAI，并測量內(nèi)中膜厚度IMT(Intima Media Thickness)。

本發(fā)明的特點及有益效果是：

水平集算法的本質(zhì)是用高一維的曲面的零等高線作為目標(biāo)邊界去分割低一維的目標(biāo)。如果從高維度看低緯度，低緯度的拓?fù)渥兓诟呔暥戎械谋憩F(xiàn)為曲面的形態(tài)變化，不會造成曲面的拓?fù)浣Y(jié)構(gòu)變化。所以水平集算法具備相當(dāng)強的低緯度的拓?fù)淇勺冃浴１景l(fā)明引入水平集方法作為ROI的初始分割算法，方法穩(wěn)定，有助于可靠結(jié)果的獲得。Markov模型可以有效的刻畫圖像中鄰域像素將的空間信息，能夠使計算機的測量結(jié)果更加準(zhǔn)確。本發(fā)明有力的支持了臨床中頸動脈超聲圖像內(nèi)中膜厚度的測量，為IMT計算機輔助測量技術(shù)的進(jìn)一步優(yōu)化發(fā)展提供了參考，對專家手動測量的方式是很好的補充。

附圖說明：

圖1算法流程圖；

圖2初始頸動脈超聲圖像；

圖3裁剪后圖像；

圖4 ROI；

圖5濾波后圖像；

圖6初始LII邊界；

圖7初始LII、MAI邊界；

圖8初始標(biāo)號場；

圖9迭代后標(biāo)號場；

圖10最終分割圖像；

圖11最終LII、MAI邊界。

具體實施方式

本發(fā)明的目的就是為了提出一種可靠、自動頸動脈超聲圖像的測量算法，使用計算機輔助測量的方式避免手動測量方式存在的缺陷?；诖四康?，針對目前的國內(nèi)外研究現(xiàn)狀，在本發(fā)明中采用水平集算法初步分割測量圖像，并引入Markov模型，將圖像分割問題轉(zhuǎn)化為每一個像素分配標(biāo)號的過程。

本發(fā)明結(jié)合超聲頸動脈圖像的特點，在所給圖像庫的測試過程中，該發(fā)明能有效的完成分割頸動脈圖像并測量IMT的工作，具有較好的理論和使用價值。

為了實現(xiàn)上述目的，本發(fā)明采用如下的技術(shù)方案：

1)圖像裁剪。初始超聲圖像的周圍分布有病人和超聲儀器的相關(guān)信息，這些信息的存在會影響后續(xù)的圖像處理步驟，因此需要裁剪移除這些部分，只保留血管部分。

2)感興趣區(qū)域(Region of Interest,ROI)提取。自動提取ROI可以避免手動分割ROI的繁瑣，實現(xiàn)自動測量IMT的目的。

3)圖像濾波。超聲圖像以實時性、可重復(fù)性、非侵入性及成本低的優(yōu)點，成為頸動脈檢查的首選成像方式。然而，回波反射帶來的斑塊噪聲不僅會減低圖像的成像質(zhì)量，更增加了計算機處理的難度，通過對圖像進(jìn)行濾波，降低噪聲的影響，獲得更加可信的結(jié)果。

4)水平集分割。采用水平集分割算法追蹤血管壁的管腔-內(nèi)膜邊界((Lumen-Intima Interface,LII))。

5)初始標(biāo)號場。根據(jù)頸動脈血管壁各層近似平行的特性，在ROI中，將血管壁內(nèi)膜邊界向下平移一定像素距離，作為中膜-外膜邊界中膜-外膜邊界(Media-Adventitia Interface,MAI)的初始狀態(tài)。然后依據(jù)這兩條邊界將ROI分割為管腔、內(nèi)中膜、外膜三個部分，并將獲得的分割圖像作為Markov模型的初始標(biāo)號場。

6)再分割。引入Markov模型，根據(jù)初始標(biāo)號場和灰度場對圖像再次分割圖像。

7)后處理。對分割圖像后處理，移除錯分小區(qū)域，獲得最終分割圖像和LII、MAI，并測量IMT。

下面結(jié)合附圖與實施例對本發(fā)明作進(jìn)一步說明。

1)圖像裁剪。在保證遠(yuǎn)端血管壁存在的條件下，截取超聲圖像中下方320×300像素大小的區(qū)域。

2)感興趣區(qū)域(Region of Interest,ROI)提取。

自動選取閾值并對其二值化。理想情況下，超聲圖像中血管的管腔和血管壁的灰度分別呈現(xiàn)出暗、亮的特點，因此根據(jù)圖像灰度的統(tǒng)計規(guī)律，采用最大類間差法自適應(yīng)獲取灰度閾值，二值化圖像。

尋找遠(yuǎn)端血管壁的近似位置。超聲圖像中血管位置的不同、斑塊噪聲的存在，使二值化的圖像可能存在多個白色區(qū)域。為了尋找與遠(yuǎn)端血管壁對應(yīng)的白色區(qū)域，對二值圖像形態(tài)學(xué)閉運算，填補圖像中的孔洞，并根據(jù)遠(yuǎn)端血管壁位于超聲圖像下方的特點，利用連通域面積較大以及質(zhì)心縱坐標(biāo)較大的原則，選擇期望的白色區(qū)域，其與黑色區(qū)域相交的上邊界可認(rèn)為是LII的近似位置。

提取ROI。正常頸動脈血管的IMT約在0.5mm到1mm的范圍內(nèi)，對應(yīng)超聲圖像中約8-16個像素點?？紤]到血管可能存在彎曲和病變的情況，我們以LII的近似位置為標(biāo)準(zhǔn)，向上取40、向下取20像素尺度，分別作為ROI的上、下邊界。在這一步驟中可以獲得兩個分別對應(yīng)原始圖像ROI。

3)圖像濾波。

本發(fā)明選取雙邊濾波算法分別對兩個ROI圖像濾波。雙邊濾波算法是一種非線性濾波方法，計算量小，它在醫(yī)學(xué)圖像的去噪處理中有較好的表現(xiàn)，并且在去除噪聲的同時能夠較好地保持圖像的邊緣信息。該算法同時考慮了圖像像素間的距離相似性和灰度值相似性，并且可以通過高斯分布來刻畫。對圖像f(r)，基于像素間距離相似性和灰度相似性的濾波輸出為：

其中g(shù)(r)是濾波后圖像，是像素ξ和r的相似性函數(shù)，為歸一化函數(shù)。σ₁是定義域方差，σ₂是值域方差；

4)水平集分割。

Step1：初始化，假設(shè)t_i時刻的水平集函數(shù)φ(m,n,t_i)為已知，m、n是平面曲線的坐標(biāo)；

Step2：求解水平集方程：獲得t_i+1時刻的水平集函數(shù)φ(m,n,t_i+1)在整個圖像區(qū)域的值，此時記為Γ_i+1的演化曲線就是φ(m,n,t_i+1)的零等值曲線，即Γ_i+1＝{φ(m,n,t_i+1)＝0}，此時的φ(m,n,t_i+1)不再是符號距離函數(shù)；

Step3：重新初始化，使用φ(m,n,t_i+1)代替下式中的初始水平集函數(shù)φ₀，迭代求解至穩(wěn)定解，仍記為φ(m,n,t_i+1)：

其中φ_ti+1為t_i+1時刻的水平集函數(shù)。

Step4：結(jié)合φ(m,n,t_i+1)的值，求解物理量的控制方程，得到t_i+1時刻的物理量的值；

Step5：重復(fù)步驟，進(jìn)入下一時刻的計算，直至停止。

5)初始標(biāo)號場。

建立初始標(biāo)號場。根據(jù)ROI圖像垂直方向邊緣圖，將初始LII向下平移一定距離(水平方向上平均梯度值最大位置)，作為MAI的初始標(biāo)記。ROI圖像被LII和MAI劃分為三部分，由上至下依次是：管腔(標(biāo)記為1)、內(nèi)中膜(標(biāo)記為2)和外膜(標(biāo)記為3)。

6)再分割。

圖像分割是基于像素的特征屬性和區(qū)域?qū)傩越o每一個像素分配標(biāo)號的過程。在馬爾可夫隨機場模型(Markov Random Field,MRF)中，常用兩個隨機場來描述待分割的圖像，一個是標(biāo)號場X，常稱為隱隨機場，用先驗分布P(X)描述標(biāo)號場的局部相關(guān)性。另一個是灰度場或特征場Y，常以標(biāo)號場為條件，用分布函數(shù)P(Y|X)描述觀測數(shù)據(jù)或特征向量的分布。其中標(biāo)號場X是馬爾可夫場，滿足Markov特性，其先驗分布可表示為

其中A＝{1,…a}是圖像中像素的集合，a是圖中像素總數(shù)，x_i∈{1,…,q}是像素i的類別標(biāo)號，q是圖像中像素類別的總數(shù)，N_i為像素i的鄰域。根據(jù)Hammersley-Clifford等價定理，所有像素的聯(lián)合概率分布為

其中Z是歸一化常數(shù)，U(x)是先驗?zāi)芰?/p>

c是一個基團，即包含若干位置的集合，x_i、x_j是像素點i、j的類別標(biāo)號，極端情況下可以認(rèn)為鄰域結(jié)構(gòu)的所有子集都是基團，C表示基團的集合。V_c(x)是定義在基團c上的勢函數(shù)，它僅與像素的鄰域像素有關(guān)，例式(10)中V₁(x_i)對應(yīng)于鄰域相關(guān)函數(shù)，勢團模型選擇二階的Potts模型進(jìn)行估計，即僅考慮二階鄰域系統(tǒng)中兩點間的相互作用。則用偽似然近似逼近的先驗概率表達(dá)為：

Potts模型：

其中β是控制鄰域間作用強度的函數(shù)，N_i、n_i是像素i的鄰域像素。

另一個隨機場Y指待分割圖像，它是可觀測的隨機過程，即P(Y)是固定的。則根據(jù)貝葉斯定理：

P(X|Y)∝P(X)P(Y|X) (13)

將模型轉(zhuǎn)化為對似然函數(shù)P(Y|X)和先驗概率P(X)的求解。這里P(X|Y)是圖像的后驗概率，其中條件分布函數(shù)P(Y|X)，通常用高斯分布描述：

其中為相同標(biāo)記(標(biāo)記為x_i)的區(qū)域像素灰度的方差，為像素i鄰域灰度的均值。

引入Markov模型，根據(jù)公式(11)、公式(14)求解模型。選取Potts模型作為公式(10)中的勢團表達(dá)函數(shù)。本發(fā)明中采用條件迭代算法(Iterative Conditional Mode,ICM)求解Markov模型。ICM是一種確定性算法，能夠獲得局部最優(yōu)解。算法在每一步的迭代過程中，都要求達(dá)到最大化以便求解過程中能過快速地收斂達(dá)到局部最優(yōu)解。ICM的計算量大大減少。

7)后處理。

我們采用最大連通域準(zhǔn)則進(jìn)行處理，移除離散小區(qū)域，得到最終分割圖像。用Sobel算子提取分割圖像中第二部分內(nèi)中膜(標(biāo)記為2)的邊界，獲取LII和MAI的最終輪廓。并計算邊界間的平均距離，作為計算機輔助的頸動脈IMT距離。