本發(fā)明屬于醫(yī)學圖像特殊組織分析與提取方法領(lǐng)域，涉及一種三維血管中軸線提取方法。
背景技術(shù)：
：目前對血管相關(guān)的醫(yī)學圖像處理，首先運用成熟的方法對相關(guān)進行分割，針對不同數(shù)據(jù)的特點運用不同的方法，針對血管造影成像進行分析處理，血管造影數(shù)據(jù)在圖像上的血管部分，往往表現(xiàn)具有“高亮”的特征，即其值相對局部范圍比較高。針對這一特點，目前對醫(yī)學圖像的三維血管中軸線提取的技術(shù)主要分為以下三大類：第一類、形態(tài)學處理方法即用形態(tài)學的腐蝕方法，對預處理過的血管數(shù)據(jù)進行三維細化。由Palágyi,Kálmán、Kuba,Attila等人提出最初的三維細化算法，而后許多研究都在此框架上提出自己的改進方法，針對不同的細化對象進行了改進。第二類、距離變換方法對體數(shù)據(jù)進行距離變換獲得距離圖，從而提取血管的中軸線。HeTS等人針對虛擬內(nèi)窺鏡的特點提出了ReliablePath的中心路徑提取算法、WanM等人和BitterI等人，提出了各自距離變換的算法，應(yīng)用在醫(yī)學圖像的拓撲細化中。第三類、模型演化法利用軸點檢測，或者上述兩種方法獲得初始的中軸線，再根據(jù)現(xiàn)有模型進行迭代，使其能量達到最優(yōu)，從而得到最合適的中軸線。K.Krissian等分析了圓柱模型的特點提出了響應(yīng)的血管分割算法并提取出中心線。此外還有大量相關(guān)研究在此類方法上有深入地探討。以上三種方法各有各的特點，形態(tài)學的方法能夠得到比較連續(xù)的中軸線，但它的結(jié)果是基于原始圖像的分割結(jié)果，若是分割結(jié)果不精準或者存在一些噪點偽影等干擾，容易生成很多不正確的或多余的軸線；距離變換方法受原圖干擾較小且具有一定的旋轉(zhuǎn)不變性，縮放不變性和平移不變性，獲得的中軸線位置準確性高，但是中軸線大多不連續(xù)，需要進行后續(xù)的處理，且在三維情況中距離變換的方法相對比二維的情況更復雜；模型法需要一個初始軸線，所以可能引入誤差，且模型的參數(shù)及演化的方法都需進一步考量，算法的復雜度高。隨著計算機技術(shù)和醫(yī)學成像技術(shù)日新月異的巨大發(fā)展，醫(yī)學圖像在對輔助醫(yī)生對病人病情診斷和治療的發(fā)揮其越來越大的作用。在醫(yī)學圖像中存在豐富的信息反映病人的身體情況，然而目前醫(yī)生只通過自身的經(jīng)驗和軟件簡單的顯示來觀察，容易受到圖像自身噪聲或者醫(yī)生主觀的判斷的影響，可能產(chǎn)生不準確診斷結(jié)果。計算機醫(yī)學圖像處理技術(shù)，作為在計算機和醫(yī)學之間的一個交叉學科，通過對醫(yī)學和成像知識和計算機圖像處理兩個方面對醫(yī)學圖像進行分析，能夠有效地利用醫(yī)學圖像反映的信息，把它以更直觀更全面的方式展現(xiàn)處理，從而達到輔助醫(yī)生進行更有效更準確地診斷的目的。由于成像方式的不同，所得到的醫(yī)學圖像數(shù)據(jù)也有各自的特點，目前研究較多的，一種是核磁共振成像MRI，它的數(shù)據(jù)分辨率高，圖像質(zhì)量較好，并且根據(jù)不同的信號加權(quán)可以獲得很多類型的成像，便于醫(yī)生和學者進行分析研究；另一種則是計算機斷層掃描成像CT，它的成像速度相比MRI要快很多，而且成本較低，各個組織成分的CT值在CT數(shù)據(jù)中都有一定的范圍，但是，它也有成像質(zhì)量一般，容易受到干擾產(chǎn)生偽影等一系列的缺點。針對不同成像方式得到的醫(yī)學圖像數(shù)據(jù)，國內(nèi)外學者用不同的方法去分析他們，達到的想要的結(jié)果。技術(shù)實現(xiàn)要素：為了克服已有三維血管中軸線提方式的準確性較差的不足，本發(fā)明提供了一種準確性較好的基于管狀特征增強濾波和脊線跟蹤的三維血管中軸線提取方法。本發(fā)明解決其技術(shù)問題所采用的技術(shù)方案是：一種基于管狀特征增強濾波和脊線跟蹤的三維血管中軸線提取方法，所述提取方法包括以下步驟：1)血管造影圖像的預處理，過程如下：1.1)最大密度投影處理四維腦部血管造影數(shù)據(jù)；1.2)多尺度管狀增強濾波；1.3)誤增強部分剔除；2)根據(jù)血管方向信息，加上血管中心值是局部最大的特點，采用類似脊線跟蹤的方法來提取三維腦補血管中軸線，過程如下：2.1)確定起始跟蹤位置：在初始化數(shù)據(jù)后選擇全局最大值做為血管跟蹤的起始位置；2.2)消除非軸點：跟蹤一根血管后，利用三維區(qū)域生長的方法，將該血管所處于的局部區(qū)域內(nèi)的非軸點體素進行消除；2.3)繼續(xù)跟蹤或者停止跟蹤：在此找到下一全局最大值做為下一根血管的初始跟蹤位置，若全局最大值低于設(shè)定的門限值，則表示整個數(shù)據(jù)已經(jīng)處理完畢，停止跟蹤。進一步，所述步驟2.1)的過程如下：2.1.1)正方向的跟蹤：根據(jù)當前點位置想pos(x0,y0,z0)和當前方向E1(a,b,c)，構(gòu)建其法平面方程a*(x-x0)+b*(y-y0)+c*(z-z0)＝0，設(shè)定一個局部范圍找到當前位置所處的局部平面；2.1.2)判斷是否處于當前平面的局部最大值，若不是，則進行局部最大偏移校正；若是，則根據(jù)血管方向E1確定下一跟蹤位置搜索區(qū)域；2.1.3)在下一跟蹤位置搜索區(qū)域中找到最大值，判斷是否是圖像邊緣或者與當前位置的值相差過大或者方向差異在一定范圍內(nèi)，若符合要求，該位置是跟蹤的下一位置nextpos則返回2.1.1)步驟繼續(xù)跟蹤；若不符合，則換另外以位置繼續(xù)判斷；如果所有位置都不符合，停止正方向的跟蹤；2.1.4)逆方向跟蹤：與正方向相似，只把跟蹤的軸向旋轉(zhuǎn)180°。再進一步，所述提取方法還包括以下步驟：3)血管中軸的分叉和連結(jié)處理：整個血管連結(jié)操作按中軸線的首部、尾部分為以下兩端點的兩種操作；對于普通連結(jié)情況：3.1)首部連結(jié)尾部的情況，2個連接點構(gòu)成新的血管，將2個中軸線合并；3.2)首部連結(jié)首部的情況，2個連接點構(gòu)成新的血管，將2個中軸線合并，且將被連結(jié)的中軸線的所有軸點方向進行反置；3.3)尾部連結(jié)首部，2個連接點構(gòu)成新的血管，將2個中軸線合并；3.4)尾部連結(jié)尾部，2個連接點構(gòu)成新的血管，將2個中軸線合并，且將被連結(jié)的中軸線的所有軸點方向進行反置；對于分叉處理情況：需補齊相的軸點；當我們找到2個端點的對應(yīng)關(guān)系時，獲得當前位置的一個方向，根據(jù)方向進行逐步跟蹤來補全中軸線。一直遍歷之前得到的初始中軸線，直到不會產(chǎn)生新的軸線的連結(jié)操作，完成全局血管中軸提取。更進一步，所述步驟1.1)中，從一個時間序列中獲得的不同時刻的造影數(shù)據(jù)，由于動脈血流在造影劑的作用下能夠在成像中反應(yīng)出高值的特性，利用最大密度投MIP進行時間維度的合并，同時分析各個時間點每個體素的值的變化率，剔除由于成像時造成的偽影的造影數(shù)據(jù)，從而獲得完整的一個血管造影數(shù)據(jù)。所述步驟1.2)中，針對不同半徑大小的血管，利用尺度空間理論，不同半徑的血管在符合它尺度下的濾波能夠取得相應(yīng)的最大響應(yīng)值，通過不同尺度的高斯濾波，構(gòu)造多尺度的管狀增強濾波。所述步驟1.3)中，通過對每一切層的梯度幅值圖像，進行最大類間分割得到梯度幅值的門限值，并且根據(jù)原數(shù)據(jù)建議的窗寬窗位值設(shè)定原數(shù)據(jù)的一個門限值，將獲得的目標區(qū)域進行剔除，目標區(qū)域即誤增強的區(qū)域，再最終結(jié)果歸一化并與原數(shù)據(jù)點乘，得到待提取中軸的數(shù)據(jù)。本發(fā)明針對人腦進行計算機斷層掃描得到的血管造影CT，提取血管中軸線并跟蹤血管走向，從而分析整個腦部血管樹，便于后續(xù)的腦部血管三維重建，血流分析等提供可靠且有效的信息。此外，血管的中軸線能夠?qū)ρ苓M一步進行精確分割或者不同模態(tài)圖像的精確配準提供了一定的拓撲先驗，能夠是這方面的工作更加有效的進行。本發(fā)明的有益效果主要表現(xiàn)在：準確性較好。附圖說明圖1是在時間序列上同一位置的橫斷面腦補血管造影原數(shù)據(jù)，最大密度投影法的運用的示意圖。圖2是不同時刻同一位置的像素變化情況的示意圖。圖3是圖2中data1、data2、data3、data4的示意圖，其中，(a)表示data1，(b)表示data2，(c)表示data3，(d)表示data4。圖4是圓柱血管模型中心體素的示意圖，E1表示最小幅值特征所對應(yīng)的方向(血管軸向)，E2和E3是剩下兩個特征向量，構(gòu)成了血管的垂直平面。圖5是不同尺度下同一橫斷面的數(shù)據(jù)，(a)為MIP的結(jié)果，(b)多尺度管狀增強濾波后進行尺度融合的結(jié)果，可以看到誤增強了最外一圈，(c)是最大尺度響應(yīng)圖，(d)是進行誤增強剔除后的結(jié)果，(e)是不同尺度下管狀增強濾波后的結(jié)果，(f)是不同尺度下的梯度幅值。圖6是血管三維跟蹤提取中軸線流程圖。圖7是血管的兩種連結(jié)情況，(a),(c)表示普通斷裂的連結(jié)，(b),(d)表示分叉的連結(jié)。具體實施方式下面結(jié)合附圖對本發(fā)明作進一步描述。參照圖1～圖7，一種基于管狀特征增強濾波和脊線跟蹤的三維血管中軸線提取方法，所述提取方法包括以下步驟：1)血管造影圖像的預處理，過程如下：1.1)最大密度投影處理四維腦部血管造影數(shù)據(jù)；1.2)多尺度管狀增強濾波；1.3)誤增強部分剔除；步驟1.1)中，從一個時間序列中獲得的不同時刻的造影數(shù)據(jù)，由于動脈血流在造影劑的作用下能夠在成像中反應(yīng)出高值的特性，利用最大密度投影(MIP)進行時間維度的合并(圖1)，同時分析各個時間點每個體素的值的變化率，剔除由于成像時造成的偽影的造影數(shù)據(jù)(圖2)，從而獲得完整的一個血管造影數(shù)據(jù)。從圖2的曲線分析圖中可以觀察到，造影部位，普通部位和偽影部位的值隨時間變化的情況，造影的部位(Data2)出現(xiàn)類似拋物線的均勻變化，普通組織部位(Data3)和背景部分(Data4)的值基本不變，偽影部位(Data1)在某一時刻產(chǎn)生突變，而且那些時刻造影劑的部分已經(jīng)失效，那些部位的CT值和普通部位接近。因此我們可以提出出現(xiàn)偽影的那些時間刻的成像，在進行最大密度投影操作。步驟1.2)中，除去原數(shù)據(jù)中非血管的部分，增強原數(shù)據(jù)中的血管部分把圖像的每個體素看成一個函數(shù)I(x,y,z)，將它進行泰勒二階展開得到I(χ0+Δχ)≈I(χ0)+ΔχT▽I(χ0)+ΔχTH(χ0)Δχ其中▽I(χ0)代表梯度，H(χ0)代表曲率，反應(yīng)了梯度的變化，每個體素經(jīng)過不同維度的二次微分都能夠得到6種二階導，它們構(gòu)成了體素的Hessian矩陣通過求解該矩陣的特征方程，能夠得到3組個特征值和它們相對應(yīng)的的特征向量，將它們按照幅值從小到大排列，表示為λ1、λ2、λ3(|λ1|＜|λ2|＜|λ3|)，通過分析這3個特征值的大小，可以找到一下規(guī)律，如表1：λ1λ2λ3結(jié)構(gòu)LLL噪聲或無結(jié)構(gòu)LLH+黑暗的片狀結(jié)構(gòu)LLH-明亮的片狀結(jié)構(gòu)LH+H+黑暗的管狀結(jié)構(gòu)LH-H-明亮的管狀結(jié)構(gòu)H+H+H+黑暗的團狀結(jié)構(gòu)H-H-H-明亮的團狀結(jié)構(gòu)表1在血管造影CT中，血管結(jié)構(gòu)在圖上表現(xiàn)出明亮的管狀結(jié)構(gòu)，我們可以利用|λ1|＜＜|λ2|＜|λ3|和λ2＜0、λ3＜0這兩個條件判讀體素是否是明亮的管狀結(jié)構(gòu)。根據(jù)Krissian等人的研究分析，絕對值最小的特征值所λ1對應(yīng)的特征向量e1代表了管狀結(jié)構(gòu)中軸線的切線方向，而其他兩個絕對值較大的特征值λ2、λ3構(gòu)成了與中軸線相垂直的平面，而且特征值的絕對值越大，代表越接近管的中心。如圖4所示。根據(jù)以上規(guī)律，F(xiàn)rangi等人研究設(shè)計了血管增強濾波，公式如下：α和β代表兩個項的敏感系數(shù)，一般取值0.5，C代表噪聲控制系數(shù)，取全局最大值的四分之一此外，針對不同半徑大小的血管，利用尺度空間理論，不同半徑的血管在符合它尺度下的濾波能夠取得相應(yīng)的最大響應(yīng)值，我們加入了多尺度的概念，通過不同尺度的高斯濾波，構(gòu)造了多尺度的管狀增強濾波經(jīng)過1.2)的操作，大部分的非血管結(jié)構(gòu)被剔除了，但是在圖像邊緣由于梯度劇烈變化的原因會造成相應(yīng)部分的誤增強，我們通過對每一切層的梯度幅值圖像，進行最大類間分割得到梯度幅值的門限值，并且根據(jù)原數(shù)據(jù)建議的窗寬窗位值設(shè)定原數(shù)據(jù)的一個門限值，將獲得的目標區(qū)域(即誤增強的區(qū)域)進行剔除，再最終結(jié)果歸一化并與原數(shù)據(jù)點乘，得到待提取中軸的數(shù)據(jù)。2)根據(jù)圖4E提供的血管方向信息，加上血管中心值是局部最大的特點，本發(fā)明設(shè)定了一種類似脊線跟蹤的方法來提取三維腦補血管中軸線，主要分為以下幾步：2.1)確定起始跟蹤位置：在初始化數(shù)據(jù)后選擇全局最大值做為血管跟蹤的起始位置；2.1.1)正方向的跟蹤：根據(jù)當前點位置想pos(x0,y0,z0)和當前方向E1(a,b,c)，構(gòu)建其法平面方程a*(x-x0)+b*(y-y0)+c*(z-z0)＝0，設(shè)定一個局部范圍找到當前位置所處的局部平面2.1.2)判斷是否處于當前平面的局部最大值，若不是，則進行局部最大偏移校正；若是，則根據(jù)血管方向E1確定下一跟蹤位置搜索區(qū)域2.1.3)在下一跟蹤位置搜索區(qū)域中找到最大值，判斷是否是圖像邊緣或者與當前位置的值相差過大或者方向差異在一定范圍內(nèi)，若符合要求，該位置是跟蹤的下一位置nextpos則返回2.1.1)步驟繼續(xù)跟蹤；若不符合，則換另外以位置繼續(xù)判斷；如果所有位置都不符合，停止正方向的跟蹤2.1.4)逆方向跟蹤：與正方向相似，只把跟蹤的軸向旋轉(zhuǎn)180°；2.2)消除非軸點：跟蹤一根血管后，利用三維區(qū)域生長的方法，將該血管所處于的局部區(qū)域內(nèi)的非軸點體素進行消除；2.3)繼續(xù)跟蹤或者停止跟蹤：在此找到下一全局最大值做為下一根血管的初始跟蹤位置，若全局最大值低于設(shè)定的門限值，則表示整個數(shù)據(jù)已經(jīng)處理完畢，停止跟蹤；流程圖如圖6。進一步，所述提取方法還包括以下步驟：3)血管中軸的分叉和連結(jié)處理：步驟2)能夠快速獲得初始的中軸線，但是在分叉處由于多尺度數(shù)據(jù)和三維區(qū)域生長終止的原因造成了斷裂，還有其他一些由于細小血管或者曲率很大的血管造成了跟蹤失效導致中軸線斷裂，針對上述情況進行分叉和連結(jié)處理，如圖6所示。整個血管連結(jié)操作按中軸線的首部、尾部分為以下兩端點的兩種操作；對于普通連結(jié)情況：3.1)首部連結(jié)尾部的情況(a)，2個連接點構(gòu)成新的血管，將2個中軸線合并；3.2)首部連結(jié)首部的情況(c)，2個連接點構(gòu)成新的血管，將2個中軸線合并，且將被連結(jié)的中軸線的所有軸點方向進行反置；3.3)尾部連結(jié)首部，2個連接點構(gòu)成新的血管，將2個中軸線合并；3.4)尾部連結(jié)尾部，2個連接點構(gòu)成新的血管，將2個中軸線合并，且將被連結(jié)的中軸線的所有軸點方向進行反置；對于分叉處理情況：無論是首部連結(jié)分叉處(另一中軸線的中間部分)還是尾部連結(jié)分叉處，都只需補齊相的軸點；當我們找到2個端點的對應(yīng)關(guān)系(即圖中五角星點和虛線箭頭所示的待連結(jié)點)時，我們可以獲得當前位置的一個方向(兩位置坐標相減)，根據(jù)方向進行逐步跟蹤來補全中軸線。整個過程操作一直遍歷之前得到的初始中軸線，直到不會產(chǎn)生新的軸線的連結(jié)操作，全局血管中軸提取就完成了。當前第1頁1 2 3