自適應增強數字減影血管造影圖像的方法

專利名稱：自適應增強數字減影血管造影圖像的方法
技術領域：
本發明屬于醫學成像技術領域，具體涉及一種自適應增強數字減影血管造影圖像的方法。
背景技術：
數字減影血管造影DSA(Digital Subtraction Angiography)技術在臨床已應用20多年，是心腦血管疾病無創診斷與介入治療手術導航的重要依據。DSA圖像處理中的一個關鍵任務就是進行圖像增強，以使血管的生理特征能夠更清楚的顯示出來，方便結構分析、運動分析、三維可視化等后續操作，以及圖像引導手術、腫瘤放射治療、治療評估等應用研究。由于X射線經過的組織厚度以及血液中的造影劑濃度不均勻，經過減影后的圖像不能完全消除人體組織引起的噪聲信號，背景和目標混雜在一起。另外，人體的組織結構和形狀很復雜，而且人與人之間有相當大的差異。因此血管的增強是一項困難的任務。為了增強血管，抑制背景，目前普遍采用各向異性平滑對圖像進行增強。Perona和Malik(參見P.Perona and J.Malik‘Scale-space and edge detectionusing anisotropic diffusion’.IEEE Trans.Pattern Anal.MachineIntell.，12(7)，629-639(1990))引入基本的各向異性擴散方程為圖像進行平滑處理。在圖像噪聲不大，并且選擇較好參數的情況下，利用各向異性平滑能夠取得不錯的效果。但是當噪聲較大時，該方法存在噪聲處理缺陷。因此Catte(參見F.Catte，P.L.Lions，J.M.Morel，and T.Coll‘Imageselective smoothing and edge detection by nonlinear diffusion’.SIAMJ.Numer.Anal.，29，182-193(1992))對其引入高斯平滑處理，使得各向異性平滑在應用時減少了噪聲的影響。高斯平滑處理的采用解決了噪聲的問題，但同時也引入了新的問題(1)在平滑噪聲的同時平滑了邊緣。(2)每一步迭代都要進行一次高斯平滑，計算速度非常緩慢。(3)隨著迭代的進行，噪聲越來越少，在每一步迭代都需要重新估計高斯函數的方差。
另外傳統的各向異性增強方法沒有確定如何選取傳導參數，大部分情況下都是通過實驗來確定。通常情況下，傳導參數過大會導致在不想執行平滑的地方產生平滑，傳導參數過小則易導致在需要平滑的區域產生增強，而且過小的傳導參數由于擴散作用不大而導致需要很多的迭代次數。不同的圖像要想找到最佳傳導參數，并不是一件容易的事情，因此有必要使得傳導參數能夠自適應地變化。因此尋找一種自適應的血管圖像增強方法一直是人們研究的一個重要方向。

發明內容
本發明的目的在于提供一種在數字減影血管造影圖像中增強血管數據的方法，該方法在整個迭代過程中避免了人為選擇傳導參數，即傳導參數能夠自適應地變化，從而能夠自適應的對DSA圖像進行增強，并降低了噪聲的影響，同時增強了邊緣，提高了圖像增強的效果。
本發明提出的在數字減影血管造影圖像中增強血管數據的方法按照下述步驟進行(1)設定迭代次數；(2)將原始數字減影血管造影圖像每個像素點(i，j)用小面模型進行擬和，并計算像素點(i，j)在擬合后的一階偏導數 和二階偏導數 (3計算每個像素點(i，j)在擬合后的Hessian矩陣，用矩陣H(i，j)表示H(i,j)=I^(i,j)xxI^(i,j)xyI^(i,j)xyI^(i,j)yy]]>并計算矩陣H(i，j)的兩個特征值λ(i，j)1，λ(i，j)2；(4)針對每個像素點(i，j)選取傳導參數T(i，j)的形式T(i,j)=(λ(i,j)12+λ(i,j)22)12]]>(5)針對每個像素點(i，j)選取傳遞函數 的形式c(|▿I^|)(i,j)=11+(|▿I^(i,j)|T(i,j))2]]>其中，|▿I^|(i,j)=I^(i,j)x2+I^(i,j)y2]]>(6)依據各個像素點的傳遞函數 采用各向異性擴散模型更新各個像素點灰度，從而更新原始圖像灰度；(7)重復步驟(2)至步驟(6)直到達到設定的迭代次數。
綜合考慮計算效率和計算精度，所述步驟(2)選用5×5的核函數對原始數字減影血管造影圖像進行擬和。
本發明利用小面模型對圖像進行擬和來達到降噪的效果，并根據圖像灰度的二階導數計算Hessian矩陣并計算傳導參數，這樣在整個迭代過程中不需要人為選擇傳導參數，使得本發明自適應的對DSA圖像進行增強，降低了噪聲的影響，并增強了邊緣，提高了增強的效果。


圖1為本發明流程圖；圖2示出了本發明與傳統各向異性擴散方法分別處理弱噪聲DSA圖像結果的對比圖，其中圖2a是從原始的弱噪聲DSA圖像選取的一部分(像素為200×200)，圖2b是采用Catte方法得到的圖像增強效果圖(傳導參數T＝30，50次迭代，高斯方差取為0.3)，圖2c是采用本發明迭代50次得到的圖像增強效果圖，圖2d是采用Catte方法得到的圖像增強效果圖(傳導參數T＝30，150次迭代，高斯方差取為0.3)，圖2e是采用本發明迭代150次得到的圖像增強效果圖；圖3示出了本方法與傳統各向異性擴散方法分別處理強噪聲DSA圖像結果的對比圖，其中圖3a是從原始的強噪聲DSA圖像選取的一部分(像素為200×200)，圖3b是采用傳統各向異性擴散Perona-Malik方法得到的圖像增強效果圖(傳導參數T＝30，50次迭代)，圖3c和圖3d是采用Catte方法得到的圖像增強效果圖(傳導參數T＝30，50次迭代，高斯方差分別取為0.3和0.4)，圖3e是利用本發明迭代50次得到的圖像增強效果圖，圖3f是采用Catte方法得到的圖像增強效果圖(傳導參數T＝30，150次迭代，高斯方差取為0.3)，圖3g是采用本發明迭代150次得到的圖像增強效果圖。
具體實施例方式
本發明包括以下步驟(1)設定迭代次數，一般選擇30到50次迭代。
(2)將原始DSA圖像每個像素點用小面模型進行擬和。
將圖像看作曲面，根據小面模型，原始圖像中每一像素點的某鄰域灰度曲面可以用離散正交多項式作曲面最佳擬合。綜合考慮計算效率和計算精度，本發明選用5×5的核函數對原始圖像進行擬和。例如，在橫坐標方向設定區間R＝{-2 -1 0 1 2}，縱坐標方向設定區間C＝{-2 -1 0 1 2}，R和C為原始圖像中某一像素點(i，j)具有元素數為5的對稱鄰域R×C內的坐標序數集，令r和c分別代表區間R和C中的值，小面模型的離散正交多項式的基底可以表示為1,r,c,r2-2,rc,c2-2,r3-175r,(r2-2)c,r(c2-2),c3-175c]]>因此，利用離散正交多項式構造的雙變量立方函數f(r，c)(i，j)來估計對應鄰域內點的灰度值為f(r,c)(i,j)=K1+K2r+K3c+K4(r2-2)+K5rc+K6(c2-2)+K7(r3-175r]]>+K8(r2-2)c+K9r(c2-2)+K10(c3-175c)]]>其中K1，K2，…Km…K10。為雙變量離散多項式的系數，由下式表示Km=Σ(r,c)∈Sgm(r,c)I(r,c)(i,j)Σ(r,c)∈Sgm2(r,c)]]>S是定義在R×C二維對稱鄰域，I(r，c)(i，j)是在某一像素點(i，j)的對稱鄰域內位置(r，c)處((r，c)∈S)的圖像灰度值，{g0(r，c)，g1(r，c)，…，g10(r，c))代表二維離散正交多項式的基底，其值分別為1,r,c,r2-2,rc,c2-2,r3-175r,(r2-2)c,r(c2-2),c3-175c.]]>小面擬和后的每一個像素點的灰度值可以表示為在擬合曲面的中心點的函數值，在此例子中，即取r＝0，c＝0時，f(0，0)(i，j)＝K1-2*K4-2*K6。像素點(i，j)在擬合后的一階偏導數 和二階偏導數 可以用像素點(i，j)對稱鄰域內的灰度估計值f(r，c)(i，j)在r＝0，c＝0處的一階和二階偏導數分別表示為I^(i,j)x=∂f(i,j)∂r,I^(i,j)y=∂f(i,j)∂c]]>I^(i,j)xx=∂2f(i,j)∂r2,I^(i,j)xy=∂2f(i,j)∂r∂c,I^(i,j)yy=∂2f(i,j)∂c2]]>其中，∂f(i,j)∂r=K2-175K7,∂f(i,j)∂c=K3-175K10]]>∂2f(i,j)∂r2=2*K4,∂2f(i,j)∂r∂c=K5,∂2f(i,j)∂c22*K6]]>(3計算每個像素點(i，j)在擬合后的Hessian矩陣，用矩陣H(i，j)表示H(i,j)=I^(i,j)xxI^(i,j)xyI^(i,j)xyI^(i,j)yy]]>并計算矩陣H(i，j)的兩個特征值λ(i，j)1，λ(i，j)2；(4)針對每個像素點(i，j)選取傳導參數T(i，j)的形式T(i,j)=(λ(i,j)12+λ(i,j)22)12]]>(5)針對每個像素點(i，j)選取傳遞函數 的形式c(|▿I^|)(i,j)=11+(|▿I^(i,j)|T(i,j))2]]>其中，|▿I^|(i,j)=I^(i,j)x2+I^(i,j)y2]]>
其中， 為像素點(i，j)梯度的模。
|▿I^|(i,j)=I^(i,j)x2+I^(i,j)y2]]>(6)依據各個像素點的傳遞函數 采用各向異性擴散模型更新各個像素點灰度，從而更新原始圖像灰度。
計算方法為In+1(i,j)=In(i,j)+Δt4dn(i,j)]]>上標n代表當前迭代次數，i和j分別代表原始圖像中橫縱坐標的取值，In(i，j)代表原始圖像的像素點(i，j)在第n次迭代后的灰度值，Δt為采樣間隔(一般取0.05-0.25)，dn(i，j)的表達形式為dn(i，j)＝cn(i，j-1)[In(i，j-1)-In(i，j)]+cn(i-1，j)[In(i-1，j)-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)]其中c(i,j)=c(|▿I^|)(i,j)]]>(7)重復步驟(2)至步驟(6)直到達到設定的迭代次數。
圖2示出了本發明與傳統各向異性擴散方法分別處理弱噪聲DSA圖像結果的對比圖。
圖2a是從原始的弱噪聲DSA圖像選取的一部分(像素為200×200)，圖2b是采用用Catte算法得到的圖像增強效果圖(傳遞參數T＝30，50次迭代，高斯方差取為0.3)，圖2d是采用Catte算法得到的圖像增強效果圖(傳遞參數T＝30，150次迭代，高斯方差取為0.3)，其結果是有的模糊了圖像邊緣，有的沒有很好的抑制噪聲。
圖2c和圖2f為本發明對原始圖像分別迭代50次和150次的效果圖，結果表明本發明對噪聲抑制，圖像邊緣增強都起到非常好的作用。
圖3示出了本發明與傳統各向異性擴散方法分別處理強噪聲DSA圖像結果的對比圖。
圖3a是從原始的強噪聲DSA圖像選取的一部分(像素為200×200)，圖3b是采用傳統各向異性擴散Perona-Malik方法得到的圖像增強效果圖(傳導參數T＝30，50次迭代)，圖3c和圖3d是采用Catte方法得到的圖像增強效果圖(傳導參數T＝30，50次迭代，高斯方差分別取為0.3和0.4)，圖3f是采用Catte方法得到的圖像增強效果圖(傳導參數T＝30，150次迭代，高斯方差取為0.3)。從效果圖可以看出，這些方法有的模糊了圖像邊緣，有的對噪聲的抑制不強，都不能達到很好的圖像增強效果。
圖3e和圖3g為本發明采用5×5區間進行小面模型擬合，分別對強噪聲DSA圖像迭代50次和150次的圖像增強效果圖，結果表明其有效降低了噪聲的影響，增強了邊緣，提高了增強效果。
權利要求
1.一種自適應增強數字減影血管造影圖像的方法，其特征是該方法按照下述步驟進行(1)設定迭代次數；(2)將原始數字減影血管造影圖像每個像素點(i，j)用小面模型進行擬和，并計算像素點(i，j)在擬合后的一階偏導數 和二階偏導數 (3)計算每個像素點(i，j)在擬合后的Hessian矩陣，用矩陣H(i，j)表示H(i,j)=I^(i,j)xxI^(i,j)xyI^(i,j)xyI^(i,j)yy]]>并計算矩陣H(i，j)的兩個特征值λ(i，j)1，λ(i，j)2；(4)針對每個像素點(i，j)選取傳導參數T(i，j)的形式T(i,j)=(λ(i,j)12+λ(i,j)22)12]]>(5)針對每個像素點(i，j)選取傳遞函數 的形式c(|ΔI^|)(i,j)=11+(|▿I^(i,j)|T(i,j))2]]>其中，|▿I^|(i,j)=I^(i,j)x2+I^(i,j)y2]]>(6)依據各個像素點的傳遞函數 采用各向異性擴散模型更新各個像素點灰度，從而更新原始圖像灰度；(7)重復步驟(2)至步驟(6)直到達到設定的迭代次數。
2.根據權利要求1所述的自適應增強數字減影血管造影圖像的方法，其特征是所述步驟(2)選用5×5的核函數對原始數字減影血管造影圖像進行擬和。
全文摘要
本發明公開了一種自適應增強數字減影血管造影圖像的方法，其步驟為①設定迭代次數；②將原始數字減影血管造影圖像用小面模型進行擬和；③計算Hessian矩陣，并計算該矩陣的兩個特征值λ
文檔編號A61B6/00GK101051384SQ20071005217
公開日2007年10月10日 申請日期2007年5月14日 優先權日2007年5月14日
發明者桑農, 張天序, 王國棟, 左崢嶸, 鐘勝, 王岳環 申請人:華中科技大學