專利名稱：一種利用序列數(shù)字減影血管造影圖像分割血管數(shù)據(jù)的方法
技術領域：
本發(fā)明屬于醫(yī)學成像技術領域，具體涉及一種在數(shù)字減影血管造影圖像中分割血管數(shù)據(jù)的方法。

背景技術：
數(shù)字減影血管造影(DSADigital Subtraction Angiography)技術在臨床中已應用20多年，是心腦血管疾病無創(chuàng)診斷與介入治療手術導航的重要依據(jù)。DSA圖像處理中的一個關鍵任務就是進行圖像分割，以使血管的生理特征能夠更清楚地顯示出來。結構分析，運動分析，三維可視化等后續(xù)操作，以及圖像引導手術，腫瘤放射治療，治療評估等應用研究都是以圖像分割為基礎的。由于X射線經(jīng)過的組織厚度以及血液中的造影劑濃度不均勻，經(jīng)過減影后的圖像不能完全消除人體組織引起的噪聲信號，背景和目標混雜在一起。另外，人體的組織結構和形狀很復雜，而且人與人之間有相當大的差異。因此血管的分割是一項困難的任務。
現(xiàn)有的血管分割技術大致可以分為6大類模式識別方法、基于模型的方法、基于跟蹤的方法、人工智能方法、神經(jīng)網(wǎng)絡方法、管狀物體檢測法，其中一些類型還可以細分(見“血管提取技術和算法的綜述”，ACMComputing Surveys，36(2)81-121，2004)。而閾值分割是最常見的直接檢測區(qū)域的分割方法，其優(yōu)點是簡單，同時對于不同類的物體灰度值或其他特征值相差很大時，它能很有效的對圖像進行分割，但它對圖像中不存在明顯灰度差異或灰度值范圍有較大重疊的圖像分割問題難以得到準確的結果。另外，由于它僅僅考慮圖像的灰度信息而不考慮圖像的空間信息，所以對噪聲和灰度不均勻很敏感。
閾值分割方法分為全局閾值分割和局部閾值分割技術。早期用得較普遍的是全值閾值分割方法，但其并不適用于某些情況。以亮目標和暗背景的圖像為例，目標的某些區(qū)域可能比背景的部分區(qū)域要暗，則任意一個確定的門限無法完全將目標和背景分離開來。也就是說，現(xiàn)有的全局分割技術很難在血管和背景灰度范圍接近的局部地方取得好的分割效果，不適用于灰度直方圖呈細窄的單峰分布的DSA圖像，很難分得連續(xù)的血管段。
而另一方面，現(xiàn)有的局部閾值分割技術通常對圖像進行分塊后，對每個子塊分別計算分割門限。由于未判斷子塊中是否存在血管，因此對不含血管的背景區(qū)域也將得到一個分割結果，這顯然是不合理的。雖然直方圖雙峰法計算分割門限時，可利用直方圖的雙峰性在一定程度上對子塊中是否存在血管作出判決，但其要求子塊內(nèi)所含血管點數(shù)和背景點數(shù)相當，這是十分苛刻很難滿足的，且當子塊較小，包含的象素數(shù)較少時，直方圖雙峰性的判斷本身也存在困難。
Peng和Li(見“Knowledge-based Adaptive Thresholding Segmentation ofDigital Subtraction Angiography Images”，Image and Vision Computing25(2007)1236-1270)提出了基于繁忙度測度的局部閾值分割方法，在血管直徑先驗知識的指導下將圖像分成適當大小的若干子塊，采用有重疊的圖像分塊技術對各子塊進行局部閾值分割，再融合各子塊分割結果，充分考慮了DSA圖像的灰度連續(xù)性，因而分割效果較好但是檢測不到部分微細血管，也不能達到較好的去噪效果。


發(fā)明內(nèi)容
本發(fā)明的目的在于提供一種利用序列數(shù)字減影血管造影圖像分割血管數(shù)據(jù)的方法，該方法可以有效地提取出更完整的血管結構，并且操作簡便，工作效率高。
本發(fā)明提供的利用序列數(shù)字減影血管造影圖像分割血管數(shù)據(jù)的方法，其步驟包括 (1)選取一個圖像序列，記為I(x，y，t)，該圖像序列I(x，y，t)包括開始注入造影劑到擴散到血管中的過程中的所有圖像，并且該一系列圖像為在時間上進行配準以后的數(shù)字減影血管造影圖像； (2)計算圖像序列I(x，y，t)中的每一點在其時間域上的灰度標準差S(x，y)’S(x，y)是點(x，y)的時域統(tǒng)計方差； (3)將圖像中每一個點的灰度標準差S(x，y)映射到0～255的灰度空間，得到一幅血管樹特征圖像； (4)對步驟(3)得到的血管樹特征圖像采用拉普拉斯濾波器進行銳化處理，增強邊緣； (5)對銳化后的血管樹特征圖像根據(jù)公式 求得n次迭代后的灰度值I(i，j)，迭代次數(shù)n初始值的取值范圍為30～50，根據(jù)處理后圖像血管區(qū)域被增強，背景區(qū)域被削弱的標準，不斷調(diào)整n的取值，得到最佳的圖像處理結果，作為各向異性平滑后的血管特征圖；其中，In(i，j)是第n次迭代后像素點(i，j)處的灰度值，Δt是迭代的間隔時間，dn(i，j)采用下式計算 dn(i，j)＝cn(i，j-1)[In(i，j-1)-In(i，j)]+cn(i-1，j)[In(i-1)-In(i，j)]+cn(i，j+1)[In(i，j+1)-In(i，j)]+cn(i+1，j)[In(i+1，j)-In(i，j)] Ix是點(i，j)處的x方向的邊緣梯度值，Iy是點(i，j)處y方向的邊緣梯度值； (6)采用基于繁忙度測度的局部閾值分割方法對步驟(5)各向異性平滑后的血管特征圖進行分割，得到二值圖像； (7)對步驟(6)分割后的得到的二值圖像，采用標記連通域的方法，提取二值圖像中面積最大的連通域，得到包含完整血管樹結構的分割結果。
本發(fā)明針對現(xiàn)有全局閾值分割技術以及部分局部閾值分割技術對DSA圖像處理的缺陷，根據(jù)DSA圖像的特點，提出一種基于在時間上配準后的序列圖像，利用統(tǒng)計特性方差的分割方法。本發(fā)明根據(jù)序列圖像在時間上的變化，利用了統(tǒng)計量方差來描述圖像上每一點的像素灰度值在時間上的變化信息，通過計算每一個像素點的方差值，然后映射到0～255的灰度空間，得到一幅描述每一點像素變化的特征圖像，這樣在時間上像素灰度值變化比較大的區(qū)域，即目標血管區(qū)域，方差值會比較大，反映在特征圖上，目標被增強了；而在時間上像素灰度值變化比較小的區(qū)域，即不含血管的區(qū)域，方差值會比較小，反映在特征圖上，背景被削弱了。本發(fā)明在得到了特征圖以后，采用了基于血管直徑的先驗知識的重疊分塊的局部閾值分割技術，對特征圖進行二值分割。本發(fā)明還采用了最大連通域的標記技術，利用血管樹在結構上的連通性，將血管樹較完整得分割出來。



圖1為本方法的流程圖； 圖2為本發(fā)明實施例中某一序列圖像，象素點的時間—灰度曲線，其中，圖2a和圖2b分別代表背景區(qū)域象素點和血管區(qū)域象素點； 圖3為基于象素方差值映射后得到的特征圖； 圖4為對特征圖進行拉普拉斯銳化后的結果圖； 圖5為對圖4進行各向異性平滑后的結果圖； 圖6為在典型實施例中采用繁忙度測度的分割算法的結果圖； 圖7為對分割后的二值圖像進行連通域標記后，得到的血管樹； 圖8為只采用序列中的最后一幀圖片，進行與特征圖相同的處理，首先拉普拉斯銳化，各向異性平滑，然后再用繁忙度算法分割以及連通域標記后的對比結果圖。

具體實施例方式 下面結合附圖來說明本發(fā)明的典型實施例 本發(fā)明實例包括以下步驟 (1)選取一個圖像序列，記為I(x，y，t)，該圖像序列I(x，y，t)包括開始注入造影劑到擴散到血管中的過程中的所有圖像，并且該一系列圖像為在時間上進行配準以后的數(shù)字減影血管造影圖像； (2)計算圖像上的每一點在時間域上的灰度標準差值，計算公式為S(x，y)其中是點(x，y)的時域統(tǒng)計標準差值，T是時間序列的長度，I(x，y，t)是第t幀圖像中點(x，y)的灰度值，

是點(x，y)在整個時間序列中的灰度平均值； (3)計算出來的灰度標準差值范圍大致在0～61.7之間，通過直線映射的方式將灰度標準差值變換到0～255的灰度空間中，形成了所要得到的血管樹特征圖，如圖3所示； (4)對特征圖采用一個拉普拉斯的銳化濾波器進行處理，即用一個3×3的矩陣，與圖像上的每一個點進行卷積相乘，將模塊的中心與像素點對應，將乘積的結果賦給當前的像素點； (5)對整幅圖像采用各向異性平滑來消除噪聲，根據(jù)公式，求得n次迭代后的邛度值I(i，j)，迭代次數(shù)n初始值的取值范圍為30～50，根據(jù)處理后圖像血管區(qū)域被增強，背景區(qū)域被削弱的標準，不斷調(diào)整n的取值，得到最佳的圖像處理結果，In(i，j)是第n次迭代后像素點(i，j)處的灰度值，Δt是迭代的間隔時間(一般取Δt＝0.05)，迭代次數(shù)n一般取30～50，dn(i，j)采用下式計算 dn(i，j)＝cn(i，j-1)[In(i，j-1)-In(i，j)]+cn(i-1，j)[In(i-1)-In(i，j)]+cn(i，j+1)[In(i，j+1)-In(i，j)]+cn(i+1，j)[In(i+1，j)-In(i，j)] Ix是點(i，j)處的x方向的邊緣梯度值，Iy是點(i，j)處y方向的邊緣梯度 圖像邊緣梯度較大時，系數(shù)c(x，y)較小，擴散較弱，保持邊緣；邊緣梯度較小時，擴散系數(shù)大，擴散強度大，圖像被模糊，從而實現(xiàn)各向異性的擴散能力； (6)采用Peng和Li提出的基于繁忙度測度的局部閾值分割方法，對步驟(5)各向異性平滑后的血管特征圖進行分割，得到二值圖像。具體實施步驟如下 (6.1)將DSA圖像劃分成N個互不重疊的a1×a2大小的子塊E1、E2、……、EN，其中，d≤a1≤2d，d≤a2≤2d，d為血管直徑的最大值；b1l、b2l、b3l、b4l為大小為b1×b2且包含子塊E1在內(nèi)的四個區(qū)域，其中，max(a1，a2)≤b1≤3max(a1，a2)，max(a1，a2)≤b2≤3max(a1，a2)； (6.2)子區(qū)域b1l、b2l、b3l、b4l分別根據(jù)OTSU分割方法(最大類間房差)，求出每個區(qū)域的分割閾值Tλl，并對四個子區(qū)域進行分割； 在對區(qū)域bkl分割后所得的二值圖對應的象素強度由fkl(i，j)來表示 (6.3)采用繁忙度作為判斷子區(qū)域分割是否合理的測度，來接受或者拒絕對子區(qū)域分割的假設門限。首先，利用如下公式
然后，利用如下公式計算血管存在性標準的繁忙度門限Hσ Hσ＝k·busymax+(1-k)·busymin 其中，busymax、busymin表示理想情況下b1×b2大小的區(qū)域繁忙度的最大或最小值，k是一個取值范圍為
的系數(shù)； 如果busy(bkl)＞Hσ，則判定該區(qū)域bkl只包含背景，將二值區(qū)域各象素值fkl(i，j)置為背景所對應的邏輯值；否則，即判定該區(qū)域包含血管，保留二值區(qū)域的結果，即fkl(i，j)值保持不變； (6.4)如果四個區(qū)域bkl均被判定為背景，其中k＝1、2、3和4，則將這四個區(qū)域的重疊子塊El判為背景；如果四個區(qū)域中只有一個區(qū)域bλl被判定包含血管，其中1≤λ≤4，則將該區(qū)域的分割門限Tλl作為子塊El的分割門限；如果四個區(qū)域中不只一個區(qū)域被判定包含血管，則對于針對不同區(qū)域bkl已求出的門限Tkl，選取其中最優(yōu)的一個門限作為El的分割門限Tl； (6.5)對子塊El重復步驟(6.2)-(6.4)，完成整幅DSA圖像的分割，得到分割后的二值圖像 (7)對分割以后的圖像，采用標記連通域的方法，提取二值圖像中面積最大的連通域，得到包含較為完整血管樹的分割結果。這里我們采用的是4-連通標記。
首先對分割得到的二值圖像從左向右、從上向下進行掃描，如果當前像素值為0，就移到下一個掃描位置。如果當前像素值為1，就檢查它左邊和上邊的2個近鄰像素，如果他們?nèi)紴?，就給當前像素一個新的標記，如果上述兩個近鄰元素有一個像素值為1，就把該像素的標記賦給當前像素，如果他們的值都為1，且具有相同的標記，就將給標記賦給當前像素，如果他們的值為1但是不具有相同的標記，就將其中的1標記賦給當前像素并作個記號表明這2個標記等價。在掃描終結時將所有等價的標記對歸入等價組，對各個組賦1個不同的標記，然后二次掃描圖像，將每個標記用它所在等價組的標記代替。標記后得到的血管樹的結果圖如圖7所示。
利用本方法可以得到較為完整的血管樹的結構，完成將血管分割出來的工作，為隨后的三維重建作基礎。
通常，我們都是采用某一幅相對血管比較清晰的圖像用來分割，但是由于造影劑隨血液的流動，不同部位血管中的造影劑濃度不同，使得單幀中包含的血管樹不連續(xù)，致使分割出來的血管不完整。本發(fā)明利用了造影劑隨時間不斷在血管中流動，血管區(qū)域的像素灰度值不斷變化的信息，將完整的血管結構通過像素點的標準差信息反映出來，得到了一幅富含完整血管結構信息的特征圖。對特征圖進行拉普拉斯銳化增強邊緣，用各向異性平滑消除噪聲，然后用基于繁忙度測度的分割算法進行分割，對分割得到的二值圖像進行連通域標記，最終將血管樹分割出來。該方法相比對于單幀圖像進行相同的處理，能夠更好地分割出完整的血管樹，取得了很好的效果。
權利要求
1.一種利用序列數(shù)字減影血管造影圖像分割血管數(shù)據(jù)的方法，其步驟包括
(1)選取一個圖像序列，記為I(x，y，t)，該圖像序列I(x，y，t)包括開始注入造影劑到擴散到血管中的過程中的所有圖像，并且該一系列圖像為在時間上進行配準以后的數(shù)字減影血管造影圖像；
(2)計算圖像序列I(x，y，t)中的每一點在其時間域上的灰度標準差S(x，y)，S(x，y)是點(x，y)的時域統(tǒng)計方差；
(3)將圖像中每一個點的灰度標準差S(x，y)映射到0～255的灰度空間，得到一幅血管樹特征圖像；
(4)對步驟(3)得到的血管樹特征圖像采用拉普拉斯濾波器進行銳化處理，增強邊緣；
(5)對銳化后的血管樹特征圖像根據(jù)公式求得n次迭代后的灰度值I(i，j)，迭代次數(shù)n初始值的取值范圍為30～50，根據(jù)處理后圖像血管區(qū)域被增強，背景區(qū)域被削弱的標準，不斷調(diào)整n的取值，得到最佳的圖像處理結果，作為各向異性平滑后的血管特征圖；其中，In(i，j)是第n次迭代后像素點(i，j)處的灰度值，Δt是迭代的間隔時間，dn(i，j)采用下式計算
dn(i，j)＝cn(i，j-1)[In(i，j-1)-In(i，j)]+cn(i-1，j)[In(i-1)-In(i，j)]
+cn(i，j+1)[In(i，j+1)-In(i，j)]+cn(i+1，j)[In(i+1，j)-In(i，j)]
Ix是點(i，j)處的x方向的邊緣梯度值，Iy是點(i，j)處y方向的邊緣梯度值；
(6)采用基于繁忙度測度的局部閾值分割方法對步驟(5)各向異性平滑后的血管特征圖進行分割，得到二值圖像；
(7)對步驟(6)分割后的得到的二值圖像，采用標記連通域的方法，提取二值圖像中面積最大的連通域，得到包含完整血管樹結構的分割結果。
全文摘要
本發(fā)明公開了一種利用序列數(shù)字減影血管造影圖像分割血管數(shù)據(jù)的方法。該方法利用統(tǒng)計量標準差來描述圖像上各點的像素灰度值在時間上的變化信息，通過計算各像素點的灰度標準差值，映射到0～255的灰度空間，得到一幅描述各點像素變化的特征圖像，這樣在時間上像素灰度值變化比較大的區(qū)域，即目標血管區(qū)域，標準差值會比較大，反映在特征圖上，目標被增強了；而在時間上像素灰度值變化比較小的區(qū)域，即不含血管的區(qū)域，方差值會比較小，反映在特征圖上，背景被削弱了。本發(fā)明采用了基于繁忙度測度的局部閾值分割方法，對特征圖進行二值分割，然后運用最大連通域的標記技術，利用血管樹在結構上的連通性，將血管樹較完整得分割出來。
文檔編號A61B6/00GK101127117SQ200710053200
公開日2008年2月20日 申請日期2007年9月11日 優(yōu)先權日2007年9月11日
發(fā)明者農(nóng) 桑, 張?zhí)煨? 曹治國, 汪春芳, 婷 郭, 王國棟 申請人:華中科技大學