本發(fā)明涉及一種頸動脈超聲圖像內中膜自動分割方法，屬于醫(yī)學圖像處理領域。

背景技術：

心血管疾病是威脅人類健康的首要疾病。動脈粥樣硬化會造成動脈血管壁增厚，造成管腔狹窄，因而是心血管疾病發(fā)生的主要原因。血管壁從內到外分為外膜層、中膜層和內膜層三個相互緊貼的膜層。內中膜厚度(intima-mediathickness,imt)是指為從管腔-內膜(lumen-intima,li)邊界到中膜-外膜(media-adventitia,ma)邊界的距離，即內膜和中膜兩部分厚度相加。imt的增加是動脈粥樣硬化的一個早期臨床癥狀。大量研究證明，頸動脈imt和心腦血管疾病存在顯著相關性，并且可以作為未來心血管事件的一個強預測指標。如今，頸動脈內中膜厚度已被看作是初步判別頸動脈粥樣硬化程度和心血管病變情況的一個重要的指標。

醫(yī)學超聲成像技術可以對頸動脈進行較清晰的成像，是一種有效的頸動脈粥樣硬化普查手段。借助于超聲成像技術，可以對頸動脈內中膜進行成像，進而可以分割出內中膜的邊界，實現(xiàn)imt的測量。然而，臨床采用的由醫(yī)生手動分割來獲取內中膜邊界的方法費時、枯燥，并且會存在較大的觀察者間的差異，因此迫切需要一種能夠自動分割內中膜邊界的圖像分割方法。

目前已有一些頸動脈內中膜分割算法，能夠實現(xiàn)較為準確的內中膜分割。然而，超聲圖像受固有的斑點噪聲影響，信噪比較低，且解剖學細節(jié)不足。嚴重的斑點噪聲甚至會掩蓋內中膜層，影響了內中膜自動分割方法成功率。本發(fā)明針對這一問題，基于蟻群優(yōu)化等算法，提出了一種具有更優(yōu)魯棒性的內中膜自動分割方法。

蟻群優(yōu)化算法(antcolonyoptimization,aco)是一種受蟻群爬行特征啟發(fā)而發(fā)展起來的一種仿生型的智能優(yōu)化算法，其特點是通過仿生自適應個體的局部最優(yōu)性共同確定問題的整體最優(yōu)解。該算法具有自學習功能和解的強搜索能力，具 有并行化、魯棒性強、正反饋的優(yōu)點，現(xiàn)已被引入到圖像分割領域當中。對于內中膜分割，可以采用蟻群算法，將ma和li邊界的輪廓提取問題轉換為最優(yōu)化問題進行求解。

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

為了替代手動分割步驟，本發(fā)明提出了一種超聲圖像內中膜自動分割方法。本發(fā)明根據(jù)內中膜結構獨特的雙層平行界面結構，提出了一種雙腿蟻群優(yōu)化算法，并結合多尺度高斯核相乘方法、邊緣檢測算子和snake模型，實現(xiàn)了內中膜的自動分割。

本發(fā)明的技術方案如下：

一種超聲圖像內中膜自動分割方法，包括以下步驟：

1)使用多尺度高斯核相乘方法抑制內中膜以外的背景信息，增強內中膜邊界，得到一個內中膜邊界凸顯的邊緣圖；

2)使用基于otsu閾值法和sobel算子的初始輪廓檢測方法，根據(jù)步驟2)得到的邊緣圖，獲取li和ma邊界的部分輪廓線段；

3)在步驟2)獲取的初始輪廓的基礎上，設計一種雙腿aco算法，將初始輪廓線段連接起來，得到一條完整的輪廓；

4)使用基于snake模型的輪廓優(yōu)化算法來進一步優(yōu)化步驟3)得到的輪廓，使最終輪廓更加平滑、連續(xù)、逼近真實邊緣。

上述步驟1)使用多尺度高斯核相乘方法抑制內中膜以外的背景信息，突顯出內中膜的li和ma邊界，得到一個邊緣圖。邊緣圖定義為兩個具有不同尺度的高斯密度核的濾波器與圖像卷積結果的乘積，并且只保留值為正的部分：

其中，

是一個小尺度的二維高斯函數(shù)，σ1取值1～3。而gσ2(y)是一個大尺 度的一維高斯函數(shù)，σ2取值10～20。

步驟2)使用基于otsu閾值法和sobel算子的初始輪廓檢測方法，主要分為以下幾個步驟：

(1)基于步驟1)得到的邊緣圖，使用otsu閾值法獲取一個二值化圖像；

(2)利用水平sobel算子從二值化圖像中提取出ma界面的上邊緣和下邊緣，及l(fā)i界面的上邊緣和下邊緣；

(3)刪去有冗余的邊緣線和有缺損的邊緣線；

(4)對于保留下的部分，取上邊緣線和下邊緣線的中點，作為每對邊緣線合并后的最終單根邊緣線。

步驟3)采用了一種雙腿aco算法，將初始輪廓線段連接起成一條完整的輪廓。其特征為，在內中膜起始點放置了兩只螞蟻，令它們從血管的左端爬向右端。其中一只螞蟻被固定放置在另一只螞蟻的上方，然后二者同時爬動，因而可以本看作是一只螞蟻的兩條腿。另外，定義通過初始輪廓檢測獲取的輪廓線段為螞蟻的強制路徑，當螞蟻爬行至這部分時，將被強制沿初始輪廓線段爬行，而線段之間的空隙將由aco算法來連接。其特征包括：

(1)轉移概率方程

螞蟻從手動選取的ps(xs,ys)點開始，從左向右爬行，每次固定向右移動1個橫坐標。在t時刻，第k只螞蟻的第l(l＝1,2)只腿從像素(x,i)移動至相鄰像素(x+1,j)的過程依據(jù)下轉移概率方程：

其中τ指信息素，η指像素的密度。α和β分別確定信息素和啟發(fā)信息的相對影響力。其中，α取值為1～3，β取值為3～5。

代表第k只螞蟻的第l只腿在下一步爬行中被允許的位置集合，它對應了像素(x,i)的右鄰域。必須指出的是，鄰域像素不能同時被兩條腿共享。因此，如果兩條腿的右鄰域出現(xiàn)重合情況，重合的像素將會被從中排除。

(2)全局更新法則

當所有螞蟻完成爬行之后，信息素按照如下公式進行更替：

τ(x,i)(t+1)＝ρτ(x,i)(t)+δτ(x,i)

其中，τ(x,i)(t)是第t次迭代時像素(x,i)處的信息素數(shù)量，τ(x,i)(t+1)是下一次迭代時像素(x,i)處的信息素數(shù)量。信息素的初始值τ0取值為0.1～1。ρ是衰減常數(shù)，用于仿真信息素的揮發(fā)，取值為0.5～1。δτ(x,i)是這一次迭代中釋放的信息素數(shù)量，其計算公式為：

其中q是一個常數(shù)，m是螞蟻的數(shù)量。q＝1，m取值為10～50。c(k)是第k只螞蟻在搜索路徑過程中的消耗函數(shù)，其定義如下：

其中in(x,y)(0≤in(xi,yi)≤1)是像素(x,i)歸一化后的灰度值。d(k)是螞蟻antk爬行路徑的終點距離之前手動選擇的終點pe(xe,ye)的距離。a是一個懲罰系數(shù)，設定了誤差距離d(k)在消耗函數(shù)中的權重，取值為1～2。

上述步驟4)使用基于snake模型的輪廓優(yōu)化算法來進一步輪廓，snake模型包含平滑能量項(smoothingenergy)、邊緣能量項(boundaryenergy)和均一能量項(uniformenergy)，并通過最小化如下能量泛函來實現(xiàn)：

其中y1(x)和y2(x)分別代表li和ma交界面的輪廓，參數(shù)μ控制平滑能量項的權重，ν控制均一能量項的權重。均一能量項連接著著li和ma兩條相互獨立的輪廓，使它們保持均一的距離。μ取值為0.1～0.3，ν取值為1～2。

本發(fā)明具有以下優(yōu)點：

本發(fā)明首先根據(jù)頸動脈超聲圖像特征，設計了一種多尺度高斯核相乘方法，對圖像進行預處理，抑制內中膜以外的背景信息，增強內中膜邊界。然后，在經(jīng)過預處理的圖像基礎上，設計了一種雙腿蟻群優(yōu)化算法，并和otsu閾值法、sobel邊緣檢測算子、蟻群優(yōu)化算法、snake模型等算法結合，實現(xiàn)了頸動脈超聲圖像內中膜的自動分割?；谂R床圖像的實驗結果表明，本方法取得了準確的分割結果，其誤差小于手動分割的觀察者間誤差；與此同時，本方法在臨床數(shù)據(jù)測試取得了98.7％的成功率，具有良好的魯棒性，并能夠應對被斑點噪聲嚴重污染的圖像。

附圖說明

圖1是本發(fā)明分割頸動脈內中膜的流程圖。

圖2是初始輪廓檢測過程中每一步的結果。

圖3是蟻群優(yōu)化算法中，螞蟻上腿(upperleg)和下腿(lowerleg)的準許位置集合的典型示例。

圖4是本發(fā)明分割結果示例圖。

具體實施方式

以下通過實施例對本發(fā)明做進一步說明，以便更好地理解本發(fā)明的技術方案。步驟如下：

1.使用多尺度高斯核相乘方法得到一個內中膜邊界凸顯的邊緣圖f(x,y)。f(x,y)計算方法為，將兩個具有不同尺度的高斯密度核的濾波器與圖像卷積結果的相乘，并且只保留值為正的部分：

其中，

是一個小尺度的二維高斯函數(shù)，本實施例中σ1取值為1.5。而gσ2(y)是一個大尺度的一維高斯函數(shù)，本實施例中σ2取值為15。

2.使用基于otsu閾值法和sobel算子的初始輪廓檢測方法，獲取li和ma邊界的部分初始輪廓線段。如圖2所示，包含以下步驟：

1)基于由上一步得到的邊緣圖(a)，使用otsu閾值法獲取一個二值化圖像，如圖(b)；

2)利用水平sobel算子從二值化圖像(b)中提取出ma界面的上邊緣和下邊緣，及l(fā)i界面的上邊緣和下邊緣；如圖(c)所示，白色線表示上邊緣，灰色線表示下邊緣。

3)經(jīng)過2)步驟后，理想結果應該含有2對(4條)輪廓線，分別對應li和ma界面的上下邊緣。然而，使用sobel算子獲取的邊緣線通常都會含有多余的錯誤邊緣線或缺口。因此，沿縱向看，多于或少于2對邊緣線的部分將會被刪去，結果如(d)所示。

4)最后，對于保留下的部分，取上邊緣線和下邊緣線的中點，作為每對邊緣線合并后的最終單根邊緣線，如(e)所示。

3.在步驟2獲取的初始輪廓的基礎上，設計一種雙腿aco算法，將初始輪廓線段連接起來，得到一條完整的輪廓。其步驟包括：

1)初始化信息素矩陣。

2)將一直螞蟻放置在手動選擇的初始點ps(xs,ys)附近(上腿放置在點上方2個像素處，下腿放置在點下方2個像素處)。

3)螞蟻從ps(xs,ys)開始向右爬動，每次固定向右移動1個橫坐標，爬行的總步數(shù)等于xe-xs。

在t時刻，第k只螞蟻的第l(l＝1,2)只腿從像素(x,i)移動至相鄰像素(x+1,j)的過程依據(jù)下轉移概率方程：

其中τ指信息素，η指像素的密度。α和β分別確定信息素和啟發(fā)信息的相對影響力。本實施例中取值為α＝1,β＝4。

代表第k只螞蟻的第l只腿在下一步爬行中被允許的位置集合，它對應了像素(x,i)的右鄰域；但是，如果兩條腿的右鄰域出現(xiàn)重合情況，重合的像素將會被從中排除。這樣可以保證螞蟻的上腿始終在下腿的上方，因而可以保證最終的輪廓不會交叉或重合。圖3展示了雙腿螞蟻爬行過程中的一個典型鄰域螞蟻此刻正從橫坐標為x的位置向x+1位置爬行，圖中淺色方格為螞蟻上腿(upperleg)和下腿(lowerleg)的準許位置集合，由于深色方格同時存在于上腿和下腿的爬行準許位置集合中，因而被從二者的集合中排除。

4)計算螞蟻路徑上像素點的信息素。

5)重復步驟2-4，直到所有的螞蟻完成了他們的爬行任務。

6)更新信息素矩陣，并計算信息素的揮發(fā)量。

信息素按照如下公式進行更替：

τ(x,i)(t+1)＝ρτ(x,i)(t)+δτ(x,i)

其中，τ(x,i)(t)是第t次迭代時像素(x,i)處的信息素數(shù)量，τ(x,i)(t+1)是下一次迭代時像素(x,i)處的信息素數(shù)量。本實施例中，信息素的初始值為τ0＝0.5。ρ是衰減常數(shù)，用于仿真信息素的會發(fā)，本實施例中選取ρ＝0.6。δτ(x,i)是這一次迭代中釋放的信息素數(shù)量，其計算公式為：

其中q是一個常數(shù)，m是螞蟻的數(shù)量。本實施例中選取q＝1，m＝20。c(k)是第k只螞蟻在搜索路徑過程中的消耗函數(shù)，其定義如下：

其中in(x,y)(0≤in(xi,yi)≤1)是像素(x,i)歸一化后的灰度值。d(k)是螞蟻antk爬行路徑的終點距離之前手動選擇的終點pe(xe,ye)的距離。a是一個懲罰系數(shù)，設定了誤差距離d(k)在消耗函數(shù)中的權重，本實施例中取a＝1。

7)重復2-6步驟，進行一定數(shù)量的迭代。本實施例中迭代次數(shù)設置為10。

8)選取具有最小消耗方程的螞蟻路徑作為最優(yōu)路徑，即為li和ma界面的輪廓。

4.使用基于snake模型的輪廓優(yōu)化算法來進一步優(yōu)化輪廓，snake模型通過最小化如下能量泛函來實現(xiàn)：

其中y1(x)和y2(x)分別代表li和ma交界面的輪廓，參數(shù)μ控制平滑能量項的權重，ν控制均一能量項的權重。本實施例選取μ＝0.1，ν＝1.4。

圖4顯示了cgaco在三張不同示例圖上的分割結果。圖a是一個噪聲嚴重的圖像。由于噪聲的污染，內中膜已經(jīng)很不清晰。如圖b所示，嚴重的斑點噪聲并沒有影響到內中膜的準確分割。圖c是一個存在彎曲情況的頸動脈圖像，圖e是內中膜出現(xiàn)典型的明顯增厚的頸動脈圖像。圖d和圖f分別對應它們的分割結果。