子宮肌瘤超聲造影圖像中的血管灌注分布提供了有用的病理和生理信息,提取出血管灌注區域有助于量化評價子宮肌瘤的血供情況。血管灌注區域像素灰度變化有別于其它區域,本文因此提出一種提取肌瘤中血管灌注區域的方法。首先對圖像序列進行去噪處理,接著采用Brox光流法對相鄰兩幀圖像進行運動估計,根據得出的位移場作warp運動校正,最后根據超聲造影圖像序列中富血供區域與乏血供區域圖像信號的差異,從周圍背景中提取出血管灌注區域。實驗結果表明,本文算法能精確地提取出血管灌注區域,可達到臨床對灌注區域識別的精度,且計算量小,實現簡單。
引用本文： 單鑫, 文銀剛, 林濤, 朱新建. 一種基于超聲造影劑的血管灌注區域的提取方法. 生物醫學工程學雜志, 2015, 32(5): 983-988. doi: 10.7507/1001-5515.20150175 復制
引言
超聲造影成像是一種通過注入造影劑來動態、清晰地顯示微細血管的新型成像技術[1],該技術可以獲得組織豐富的供血及血流灌注信息。造影微泡具有較強的散射性，注入血流后可使血液回聲增強，達到良好的“血管顯影”效果[2]。從周圍組織中提取出血管灌注區域是對造影信號定量分析的基礎,這對臨床診斷組織病變以及了解腫瘤的生長狀況具有重要臨床意義。目前對測定感興趣區域(region of interest,ROI)的時間-強度曲線(time-intensity curves，TICs)的研究較多，TICs反映了造影劑注入后組織某點增強的回波信號(即圖像中像素灰度)隨時間變化的曲線。基于TICs的特征參數進行參數成像，借助超聲灌注成像來分析血管灌注區域以此來定量分析超聲造影圖像，然而由于超聲圖像噪聲強、偽影多以及組織回波信號多樣性，難以確定合適的數學模型，大大增加了超聲灌注參數成像的難度，因此該成像技術尚處于研究階段[3]。Frinking等[4]對造影劑在不同的聲輻射條件下，對靈敏度進行解相關運算并比較不同結果，在此基礎上提出了一種利用微泡破壞機制的多脈沖釋放成像(multi-pulse release imaging)方法，對血管灌注區域進行成像。該實驗系統有兩個超聲探頭：釋放破壞性脈沖信號的釋放探頭以及用于超聲成像的探頭。這種成像方法要求釋放探頭與成像探頭共同聚焦于同一成像目標區,兩個探頭很難協調掃描,因此該方法很難實現。宋延淳等[5]提出一種基于高機械指數、單成像脈沖發射導致微泡破壞的射頻回波解相關(decorrelation,DCR)檢測方法,在仿體實驗上取得較好的實驗效果。以上是基于爆破樣成像所產生的回波信號進行解相關運算，建模及分析過程復雜，也僅實現于離體灌注仿體模型中，尚未應用于臨床數據。
本文以子宮肌瘤超聲造影圖像為研究對象，創新地提出了一種直接對造影圖像序列進行動態分析，通過造影微泡與周圍組織成像的差別提取出血管灌注區域的方法。本文采用Brox光流法估計出相鄰圖像幀的運動位移并根據位移場利用warp方法進行校正，基于校正后的肌瘤造影圖像序列，根據造影微泡在血液中表現出的高亮、閃爍,區別于其他背景組織的信號特點提取出灌注區域。實驗結果表明，本文的算法應用于臨床數據時，不但提取效果好，滿足臨床要求，而且算法復雜度小，實現簡單。
1 原理與方法
提取血管灌注區域的方法流程如圖 1所示。先將子宮肌瘤超聲造影錄像解碼為一幀幀圖像，接著對每幀圖像進行高斯低通濾波，以減少斑點噪聲影響，再對預處理后的圖像序列進行運動估計并校正，得到消除呼吸運動影響后的圖像序列，最后通過從造影圖像序列上采樣得到的分割閾值將血管灌注區域提取出來。
圖1
血管灌注區域的提取策略基本流程圖
Figure1.
Basic flow chart of the strategy for extraction of vascular perfusion region
1.1 運動估計與校正
子宮肌瘤屬于盆腔腫瘤，在超聲圖像中的運動主要由呼吸引起。為了減少呼吸運動的影響，控制呼吸是最常見的方法，然而錄制造影過程的時間往往需要20 s以上，這對患者來說是很困難的。另外可以對超聲造影圖像序列進行運動校正，由于超聲圖像本身噪聲多、偽影多，加上造影劑的不斷攝入，使得相鄰造影圖像間出現明顯的灰度差異，致使基于配準的的運動校正算法失效。因此針對該特點，我們需要尋求一種抗噪聲能力強、識別組織運動部分準確性高的運動估計算法。基于此，本文結合Brox光流法估計造影圖像相鄰幀的運動場,并根據估計出的運動場采用warp技術[6]校正造影圖像。
光流法是分析運動圖像的重要方法，光流表達了圖像的變化，包含了圖像的運動信息。光流場的計算方法最初由Horn和Schunk提出[7]。本文使用的是Brox等[8]于2004年提出的一種改進的變分光流法。
Brox光流法的能量函數為：
| $ E\left( {u,v} \right) = {E_{{\rm{Data}}}} + \alpha {E_{{\rm{Smooth}}}} $ |
展開如式(2)、(3)：
| $ \begin{array}{l} {E_{{\rm{Data}}}}\left( {u,v} \right) = \\ \int_\Omega {\psi \left( {\left| {I\left( {x + w} \right) - I{{\left( x \right)}^2}} \right| + \gamma {{\left| {\nabla I\left( {x + w} \right) - \nabla I\left( x \right)} \right|}^2}} \right){\rm{d}}} x \end{array} $ |
| $ \begin{array}{l} {E_{{\rm{Smooth}}}}\left( {u,v} \right) = \\ \int_\Omega {\psi \left( {{{\left| {{\nabla _3}u} \right|}^2} + {{\left| {{\nabla _3}v} \right|}^2}} \right){\rm{d}}} x \end{array} $ |
其中令x=(x,y,t)T,w=(u,v,1)T；表示圖像的空間度,表示時間-空間梯度，Ix、 Iy、 It分別表示像素灰度沿x、y、t(時間)方向的梯度,u和v分別表示此像素點在x和y方向的位移。EData(u,v)是灰度連續假設與梯度連續假設的能量函數，γ為兩種假設之間的權重，γ越大，算法對灰度守恒的依賴越小。ESmooth(u,v)為光流場平滑假設的能量函數，該假設是為了避免出現圖像梯度消失而導致算法估計出異常值的情況，α為數據項(2)與平滑項(3)之間的權重。為了提高了算法的魯棒性，通過二次平方函數(ε為大于0的常數)加強了對孤立點的懲罰。
最后求出使得總能量函數式(1)最小值時的(u,v)即圖像運動位移場，根據(u,v)對相鄰兩幀圖像每個像素進行了warp技術的校正。
1.2 血管灌注區域的提取
超聲造影成像是利用微泡造影劑經過外周靜脈注射后進入全身血液循環，使得病灶處的回聲和血流信號增強而提高超聲診斷的分辨力、敏感性和特異性的技術[2]。造影微泡的背向散射較紅細胞強得多，而且與組織不同的是它還是聲波發射體，在入射聲強的影響下，造影微泡會發生膨脹與收縮振動，當膨脹幅度超過收縮幅度時，導致波形畸變，產生諧波效應；當超聲功率達到某一臨界點時，微泡局部受到的壓力大于它本身可承受的壓力,造影微泡會被破壞，產生強烈的背向散射，顯示出該病灶區血管容量的信息，這些都是周圍組織并不具備的特性[9-11]。判斷子宮肌瘤內血供豐富區域主要依賴于肌瘤相對于周圍正常組織的回聲情況。子宮肌瘤周圍的組織有呈強回聲的，有呈低回聲的，而肌瘤內血供豐富的區域在灌注峰值的時候表現為強回聲，單從孤立的一幀圖像上無法區分出呈強回聲的部分是正常組織還是肌瘤血供豐富的區域。當多幀圖像連續播放時，可以發現與正常組織強回聲區域在多幀圖像中強度變化不大相同，造影微泡灌注的區域呈現出高增強并且強度會發生很大的變化，這是造影微泡在入射聲強的影響下不斷發生膨脹、壓縮甚至破滅的動態過程表現。可以把超聲成像束看成是有一定厚度的平面薄板，薄板上的造影氣泡流進、流出、破滅、波動等變化在灌注峰值時達到一個動態平衡[12]，在超聲圖像上表現為高亮、閃爍的似血液樣的流動區域，明顯區別于周圍正常組織以及腫瘤內部乏血供的部分。病灶處血管灌注達到峰值后造影劑并不會立刻消退，峰值信號強度會維持一段時間，根據TICs的峰值，采集峰值時期內連續若干秒的圖像序列作為提取血管灌注區域的圖像序列。因此本文根據圖像序列相鄰幀之間灌注區域與非灌注區域灰度值變化幅度的差異設計了一種提取血管灌注區域的算法。
為了降低噪聲對單個像素點信號強度的影響，首先對每幀圖像的每個像素點進行8鄰域平均。接著計算1~a幀中相鄰兩幀圖像間像素點灰度的平均變化值，計算公式如下：
| $ T = \frac{{\sum\limits_{s = 1}^{a - 1} {\left| {I\left( s \right) - I\left( {s + 1} \right)} \right|} }}{{a - 1}} $ |
公式(4)中s為圖片序列的起始幀，a為圖片序列的總幀數，I(s)為序號為s的圖像，｜I(s)-I(s+1)｜為相鄰兩幀圖像灰度差，為減少誤差，再將所有比較的結果求平均值，最終得到a幀圖像間灰度差的平均值矩陣T,此時設定一個灰度值變化的閾值t，矩陣T中高于t的點認為是灰度變化劇烈的點，即造影微泡灌注的位置，小于t的點認為是灰度變化平緩的點，即非血管灌注區域的點，利用閾值分割的方法即可提取出有造影微泡存在的點。其中閾值t的選取至關重要，目前還沒有關于超聲造影圖像灰度值變化方面的統計研究，本文采取的是采樣平均法。通過播放超聲造影錄像在超聲科醫生的指導下初步判斷出血管灌注區域，在該區域中選取3個很小的區域作為觀察對象，用式(4)求取樣本區域的T，分別記為：T1，為m1×n1矩陣；T2，為m2×n2 矩陣；T3，為m3×n3矩陣，求取t的公式如下：
| $ t = \frac{{\left( {\sum\limits_1^{{m_1}} {\sum\limits_1^{{n_1}} {{T_1}} } } \right)\frac{1}{{{m_1} \times {n_1}}} + \left( {\sum\limits_1^{{m_2}} {\sum\limits_1^{{n_2}} {{T_2}} } } \right)\frac{1}{{{m_2} \times {n_2}}} + \left( {\sum\limits_1^{{m_3}} {\sum\limits_1^{{n_3}} {{T_3}} } } \right)\frac{1}{{{m_3} \times {n_3}}}}}{3} $ |
最終提取出T≥t的點，標記為血管灌注區域的點。
2 結果與分析
試驗的運行環境是i5-3210M @ 2.50 GHz 雙核，Windows 7，64位，MATLAB 2011b。試驗數據來源于重慶海扶醫院，子宮肌瘤患者在高強度聚焦超聲治療前需要先靜脈推注2.0 mL(10.0 mg)造影劑聲諾維，以觀察肌瘤邊界及血流灌注情況。超聲實時監控設備為意大利百盛公司生產的 My-Lab70 型 B 超，頻率1.0~8.0 MHz，超聲探頭置于聚焦超聲換能器的中央，注射造影劑前將超聲探頭調至患者矢狀位、子宮肌瘤最大切面處。有研究表明,經肘靜脈注入造影劑至子宮肌瘤開始顯影時間，平均約為15 s,肌瘤開始灌注到達到峰值強度約為12 s[13]。本文采集數據時，待注入造影劑10 s后開始錄制造影過程，錄制時間為53 s，幀率為5 fps，獲得265幀圖像。我們需要的是子宮肌瘤灌注達到峰值時期的圖像，觀看錄像，發現錄制至第30秒的時候圖像灌注達到峰值強度并持續超過10 s，取30 s到34 s共20幀圖像，都處于灌注峰值期，編號為1~20，為了避免相鄰兩幀間造影劑的變化過于微小，取編號為奇數的圖像，形成新的圖像序列，重新編號為1~10。
先對1~10幀中兩兩相鄰圖像間進行運動估計并校正。第1、2幀圖像的估計結果如圖 2所示，圖 2(b)是圖 2(a)中白色矩形框中像素的光流場,箭頭表示位移方向，線段的長短表示位移量的大小，為顯示方便，作了放大處理。白色矩形框選取在肌瘤的邊界處，可見光流場的左上部分有明顯的位移，與觀看造影錄像觀察到的位移方位一致，右下部分光流的場的模很小，說明幾乎沒有位移。這跟實際圖像的情況是一致的：左上部分是腫瘤部分，隨著呼吸會有幾個像素的位移，而右下角為貼近皮膚部分，位移很小。
圖2
子宮肌瘤超聲造影序列及矩形框中的光流場
(a)原始子宮肌瘤超聲造影圖像序列；(b)第1幀與第2幀圖像白色矩形框中像素的光流矢量圖
Figure2. Ultrasonic imaging sequences of fibroids, and optical flow field in rectangular box(a) ultrasonic imaging sequences of fibroids; (b) optical flow vector diagram in white square window marked in picture 1 and 2
為了驗證校正運動算法的性能，本文采用計算互信息(mutual information,MI)和互相關系數(cross correlation,CC)作為衡量指標[14]，結果如表 1所示。
表1
兩兩相鄰幀圖像運動校正結果
Table1.
Motion correction results of two neighboring frame images
MI是兩幅圖像相關性的測度,MI越大說明兩幅圖像越接近,CC是用來表達兩幅圖像的線性相關性,該值越大表明圖像校正的精確度越高。從表 1可以看出校正后相對于校正前，MI增大，CC增大，表明校正算法的有效性。
確定灌注區域的分割閾值t，需要先在圖像中選取3個很小矩形區域，作為采樣區域，即T1、T2、T3，如圖 3所示。
圖3
采樣區域T1、T2、T3
Figure3.
Sampling area T1, T2, and T3
接著先由式(4)處理校正后的圖像，再根據式(5)，由T1、T2、T3計算得t=7.699 6。 提取結果如圖 4所示，提取出來的點處于腫瘤內部，超聲遠場以及腫瘤周邊的點幾乎排除在灌注區域之外。
圖4
血管灌注區域提取結果
(a)二值圖像表示提取結果；(b)在造影圖像上用紅色點標記出灌注區域，藍色曲線內為ROI
Figure4. Results of vascular perfusion region extraction(a) binary image of extraction results; (b) red points: perfusion area; turquoise curve: ROI
分別從來自校正后的和對灌注區域處理后的造影圖像序列的ROI[圖 4(b)中藍色曲線中的區域]中提取TICs，用于評估提取算法的性能[15]。錄制的53 s造影錄像共265幀圖像，都兩兩經過了warp校正，其中第20 s起造影微泡開始進入肌瘤。根據先前提取出的血管灌注區域的像素點坐標，將20~53 s采集到的圖像中相同位置點的灰度值用前10 s圖像序列對應位置點的像素值代替，也就是保留灌注區域的基礎像素值，消除造影微泡對像素灰度值的影響。提取出的TICs如圖 5所示。前100幀圖像的兩組TICs是重合的，因為此時造影劑還未進入肌瘤，從100幀起造影劑開始進入肌瘤內部，對灌注區域處理后的圖像組的TICs信號強度并未發生劇烈變化，僅有輕微的增強，表明被替換的點是有造影微泡存在的點，即算法準確地識別出了血管灌注區域。
圖5
圖像序列
紅色曲線為運動校正后圖像序列的ROI的TICs，藍色曲線為對灌注區域處理后的圖像序列ROI的TICs
Figure5. Image sequencesred curve： the TICs of ROI that from after correction； blue curve： TICs of ROI from image sequences after processing infusion area
為了進一步驗證提取算法的實用性，由兩位經驗豐富的超聲科醫生觀看造影錄像，再對比由本算法提取出來的血管灌注區域，提取結果和實際灌注區域吻合度高，達到了臨床對灌注區域的識別精度，表明該算法能夠精確地提取出血管灌注區域，具有很好的實用性和準確性。
3 結論
針對子宮肌瘤超聲造影圖像的血管灌注區域的提取，本文提出了一種基于運動校正后的圖像序列閾值提取法。預處理圖像后，首先利用Brox光流法估計出運動場，再根據運動場對相鄰兩幀圖像進行warp校正，接著在腫瘤灌注區域采樣求得分割閾值t,最后根據超聲圖像中富血供與乏血供區域信號強度的變化由閾值t提取出血管灌注區域。結果表明該算法能很好地提取出腫瘤中的血管灌注區域，為臨床對造影區域的識別提供可靠的參考，為進一步超聲造影圖像的定量分析打下基礎。
引言
超聲造影成像是一種通過注入造影劑來動態、清晰地顯示微細血管的新型成像技術[1],該技術可以獲得組織豐富的供血及血流灌注信息。造影微泡具有較強的散射性，注入血流后可使血液回聲增強，達到良好的“血管顯影”效果[2]。從周圍組織中提取出血管灌注區域是對造影信號定量分析的基礎,這對臨床診斷組織病變以及了解腫瘤的生長狀況具有重要臨床意義。目前對測定感興趣區域(region of interest,ROI)的時間-強度曲線(time-intensity curves，TICs)的研究較多，TICs反映了造影劑注入后組織某點增強的回波信號(即圖像中像素灰度)隨時間變化的曲線。基于TICs的特征參數進行參數成像，借助超聲灌注成像來分析血管灌注區域以此來定量分析超聲造影圖像，然而由于超聲圖像噪聲強、偽影多以及組織回波信號多樣性，難以確定合適的數學模型，大大增加了超聲灌注參數成像的難度，因此該成像技術尚處于研究階段[3]。Frinking等[4]對造影劑在不同的聲輻射條件下，對靈敏度進行解相關運算并比較不同結果，在此基礎上提出了一種利用微泡破壞機制的多脈沖釋放成像(multi-pulse release imaging)方法，對血管灌注區域進行成像。該實驗系統有兩個超聲探頭：釋放破壞性脈沖信號的釋放探頭以及用于超聲成像的探頭。這種成像方法要求釋放探頭與成像探頭共同聚焦于同一成像目標區,兩個探頭很難協調掃描,因此該方法很難實現。宋延淳等[5]提出一種基于高機械指數、單成像脈沖發射導致微泡破壞的射頻回波解相關(decorrelation,DCR)檢測方法,在仿體實驗上取得較好的實驗效果。以上是基于爆破樣成像所產生的回波信號進行解相關運算，建模及分析過程復雜，也僅實現于離體灌注仿體模型中，尚未應用于臨床數據。
本文以子宮肌瘤超聲造影圖像為研究對象，創新地提出了一種直接對造影圖像序列進行動態分析，通過造影微泡與周圍組織成像的差別提取出血管灌注區域的方法。本文采用Brox光流法估計出相鄰圖像幀的運動位移并根據位移場利用warp方法進行校正，基于校正后的肌瘤造影圖像序列，根據造影微泡在血液中表現出的高亮、閃爍,區別于其他背景組織的信號特點提取出灌注區域。實驗結果表明，本文的算法應用于臨床數據時，不但提取效果好，滿足臨床要求，而且算法復雜度小，實現簡單。
1 原理與方法
提取血管灌注區域的方法流程如圖 1所示。先將子宮肌瘤超聲造影錄像解碼為一幀幀圖像，接著對每幀圖像進行高斯低通濾波，以減少斑點噪聲影響，再對預處理后的圖像序列進行運動估計并校正，得到消除呼吸運動影響后的圖像序列，最后通過從造影圖像序列上采樣得到的分割閾值將血管灌注區域提取出來。
圖1
血管灌注區域的提取策略基本流程圖
Figure1.
Basic flow chart of the strategy for extraction of vascular perfusion region
1.1 運動估計與校正
子宮肌瘤屬于盆腔腫瘤，在超聲圖像中的運動主要由呼吸引起。為了減少呼吸運動的影響，控制呼吸是最常見的方法，然而錄制造影過程的時間往往需要20 s以上，這對患者來說是很困難的。另外可以對超聲造影圖像序列進行運動校正，由于超聲圖像本身噪聲多、偽影多，加上造影劑的不斷攝入，使得相鄰造影圖像間出現明顯的灰度差異，致使基于配準的的運動校正算法失效。因此針對該特點，我們需要尋求一種抗噪聲能力強、識別組織運動部分準確性高的運動估計算法。基于此，本文結合Brox光流法估計造影圖像相鄰幀的運動場,并根據估計出的運動場采用warp技術[6]校正造影圖像。
光流法是分析運動圖像的重要方法，光流表達了圖像的變化，包含了圖像的運動信息。光流場的計算方法最初由Horn和Schunk提出[7]。本文使用的是Brox等[8]于2004年提出的一種改進的變分光流法。
Brox光流法的能量函數為：
| $ E\left( {u,v} \right) = {E_{{\rm{Data}}}} + \alpha {E_{{\rm{Smooth}}}} $ |
展開如式(2)、(3)：
| $ \begin{array}{l} {E_{{\rm{Data}}}}\left( {u,v} \right) = \\ \int_\Omega {\psi \left( {\left| {I\left( {x + w} \right) - I{{\left( x \right)}^2}} \right| + \gamma {{\left| {\nabla I\left( {x + w} \right) - \nabla I\left( x \right)} \right|}^2}} \right){\rm{d}}} x \end{array} $ |
| $ \begin{array}{l} {E_{{\rm{Smooth}}}}\left( {u,v} \right) = \\ \int_\Omega {\psi \left( {{{\left| {{\nabla _3}u} \right|}^2} + {{\left| {{\nabla _3}v} \right|}^2}} \right){\rm{d}}} x \end{array} $ |
其中令x=(x,y,t)T,w=(u,v,1)T；表示圖像的空間度,表示時間-空間梯度，Ix、 Iy、 It分別表示像素灰度沿x、y、t(時間)方向的梯度,u和v分別表示此像素點在x和y方向的位移。EData(u,v)是灰度連續假設與梯度連續假設的能量函數，γ為兩種假設之間的權重，γ越大，算法對灰度守恒的依賴越小。ESmooth(u,v)為光流場平滑假設的能量函數，該假設是為了避免出現圖像梯度消失而導致算法估計出異常值的情況，α為數據項(2)與平滑項(3)之間的權重。為了提高了算法的魯棒性，通過二次平方函數(ε為大于0的常數)加強了對孤立點的懲罰。
最后求出使得總能量函數式(1)最小值時的(u,v)即圖像運動位移場，根據(u,v)對相鄰兩幀圖像每個像素進行了warp技術的校正。
1.2 血管灌注區域的提取
超聲造影成像是利用微泡造影劑經過外周靜脈注射后進入全身血液循環，使得病灶處的回聲和血流信號增強而提高超聲診斷的分辨力、敏感性和特異性的技術[2]。造影微泡的背向散射較紅細胞強得多，而且與組織不同的是它還是聲波發射體，在入射聲強的影響下，造影微泡會發生膨脹與收縮振動，當膨脹幅度超過收縮幅度時，導致波形畸變，產生諧波效應；當超聲功率達到某一臨界點時，微泡局部受到的壓力大于它本身可承受的壓力,造影微泡會被破壞，產生強烈的背向散射，顯示出該病灶區血管容量的信息，這些都是周圍組織并不具備的特性[9-11]。判斷子宮肌瘤內血供豐富區域主要依賴于肌瘤相對于周圍正常組織的回聲情況。子宮肌瘤周圍的組織有呈強回聲的，有呈低回聲的，而肌瘤內血供豐富的區域在灌注峰值的時候表現為強回聲，單從孤立的一幀圖像上無法區分出呈強回聲的部分是正常組織還是肌瘤血供豐富的區域。當多幀圖像連續播放時，可以發現與正常組織強回聲區域在多幀圖像中強度變化不大相同，造影微泡灌注的區域呈現出高增強并且強度會發生很大的變化，這是造影微泡在入射聲強的影響下不斷發生膨脹、壓縮甚至破滅的動態過程表現。可以把超聲成像束看成是有一定厚度的平面薄板，薄板上的造影氣泡流進、流出、破滅、波動等變化在灌注峰值時達到一個動態平衡[12]，在超聲圖像上表現為高亮、閃爍的似血液樣的流動區域，明顯區別于周圍正常組織以及腫瘤內部乏血供的部分。病灶處血管灌注達到峰值后造影劑并不會立刻消退，峰值信號強度會維持一段時間，根據TICs的峰值，采集峰值時期內連續若干秒的圖像序列作為提取血管灌注區域的圖像序列。因此本文根據圖像序列相鄰幀之間灌注區域與非灌注區域灰度值變化幅度的差異設計了一種提取血管灌注區域的算法。
為了降低噪聲對單個像素點信號強度的影響，首先對每幀圖像的每個像素點進行8鄰域平均。接著計算1~a幀中相鄰兩幀圖像間像素點灰度的平均變化值，計算公式如下：
| $ T = \frac{{\sum\limits_{s = 1}^{a - 1} {\left| {I\left( s \right) - I\left( {s + 1} \right)} \right|} }}{{a - 1}} $ |
公式(4)中s為圖片序列的起始幀，a為圖片序列的總幀數，I(s)為序號為s的圖像，｜I(s)-I(s+1)｜為相鄰兩幀圖像灰度差，為減少誤差，再將所有比較的結果求平均值，最終得到a幀圖像間灰度差的平均值矩陣T,此時設定一個灰度值變化的閾值t，矩陣T中高于t的點認為是灰度變化劇烈的點，即造影微泡灌注的位置，小于t的點認為是灰度變化平緩的點，即非血管灌注區域的點，利用閾值分割的方法即可提取出有造影微泡存在的點。其中閾值t的選取至關重要，目前還沒有關于超聲造影圖像灰度值變化方面的統計研究，本文采取的是采樣平均法。通過播放超聲造影錄像在超聲科醫生的指導下初步判斷出血管灌注區域，在該區域中選取3個很小的區域作為觀察對象，用式(4)求取樣本區域的T，分別記為：T1，為m1×n1矩陣；T2，為m2×n2 矩陣；T3，為m3×n3矩陣，求取t的公式如下：
| $ t = \frac{{\left( {\sum\limits_1^{{m_1}} {\sum\limits_1^{{n_1}} {{T_1}} } } \right)\frac{1}{{{m_1} \times {n_1}}} + \left( {\sum\limits_1^{{m_2}} {\sum\limits_1^{{n_2}} {{T_2}} } } \right)\frac{1}{{{m_2} \times {n_2}}} + \left( {\sum\limits_1^{{m_3}} {\sum\limits_1^{{n_3}} {{T_3}} } } \right)\frac{1}{{{m_3} \times {n_3}}}}}{3} $ |
最終提取出T≥t的點，標記為血管灌注區域的點。
2 結果與分析
試驗的運行環境是i5-3210M @ 2.50 GHz 雙核，Windows 7，64位，MATLAB 2011b。試驗數據來源于重慶海扶醫院，子宮肌瘤患者在高強度聚焦超聲治療前需要先靜脈推注2.0 mL(10.0 mg)造影劑聲諾維，以觀察肌瘤邊界及血流灌注情況。超聲實時監控設備為意大利百盛公司生產的 My-Lab70 型 B 超，頻率1.0~8.0 MHz，超聲探頭置于聚焦超聲換能器的中央，注射造影劑前將超聲探頭調至患者矢狀位、子宮肌瘤最大切面處。有研究表明,經肘靜脈注入造影劑至子宮肌瘤開始顯影時間，平均約為15 s,肌瘤開始灌注到達到峰值強度約為12 s[13]。本文采集數據時，待注入造影劑10 s后開始錄制造影過程，錄制時間為53 s，幀率為5 fps，獲得265幀圖像。我們需要的是子宮肌瘤灌注達到峰值時期的圖像，觀看錄像，發現錄制至第30秒的時候圖像灌注達到峰值強度并持續超過10 s，取30 s到34 s共20幀圖像，都處于灌注峰值期，編號為1~20，為了避免相鄰兩幀間造影劑的變化過于微小，取編號為奇數的圖像，形成新的圖像序列，重新編號為1~10。
先對1~10幀中兩兩相鄰圖像間進行運動估計并校正。第1、2幀圖像的估計結果如圖 2所示，圖 2(b)是圖 2(a)中白色矩形框中像素的光流場,箭頭表示位移方向，線段的長短表示位移量的大小，為顯示方便，作了放大處理。白色矩形框選取在肌瘤的邊界處，可見光流場的左上部分有明顯的位移，與觀看造影錄像觀察到的位移方位一致，右下部分光流的場的模很小，說明幾乎沒有位移。這跟實際圖像的情況是一致的：左上部分是腫瘤部分，隨著呼吸會有幾個像素的位移，而右下角為貼近皮膚部分，位移很小。
圖2
子宮肌瘤超聲造影序列及矩形框中的光流場
(a)原始子宮肌瘤超聲造影圖像序列；(b)第1幀與第2幀圖像白色矩形框中像素的光流矢量圖
Figure2. Ultrasonic imaging sequences of fibroids, and optical flow field in rectangular box(a) ultrasonic imaging sequences of fibroids; (b) optical flow vector diagram in white square window marked in picture 1 and 2
為了驗證校正運動算法的性能，本文采用計算互信息(mutual information,MI)和互相關系數(cross correlation,CC)作為衡量指標[14]，結果如表 1所示。
表1
兩兩相鄰幀圖像運動校正結果
Table1.
Motion correction results of two neighboring frame images
MI是兩幅圖像相關性的測度,MI越大說明兩幅圖像越接近,CC是用來表達兩幅圖像的線性相關性,該值越大表明圖像校正的精確度越高。從表 1可以看出校正后相對于校正前，MI增大，CC增大，表明校正算法的有效性。
確定灌注區域的分割閾值t，需要先在圖像中選取3個很小矩形區域，作為采樣區域，即T1、T2、T3，如圖 3所示。
圖3
采樣區域T1、T2、T3
Figure3.
Sampling area T1, T2, and T3
接著先由式(4)處理校正后的圖像，再根據式(5)，由T1、T2、T3計算得t=7.699 6。 提取結果如圖 4所示，提取出來的點處于腫瘤內部，超聲遠場以及腫瘤周邊的點幾乎排除在灌注區域之外。
圖4
血管灌注區域提取結果
(a)二值圖像表示提取結果；(b)在造影圖像上用紅色點標記出灌注區域，藍色曲線內為ROI
Figure4. Results of vascular perfusion region extraction(a) binary image of extraction results; (b) red points: perfusion area; turquoise curve: ROI
分別從來自校正后的和對灌注區域處理后的造影圖像序列的ROI[圖 4(b)中藍色曲線中的區域]中提取TICs，用于評估提取算法的性能[15]。錄制的53 s造影錄像共265幀圖像，都兩兩經過了warp校正，其中第20 s起造影微泡開始進入肌瘤。根據先前提取出的血管灌注區域的像素點坐標，將20~53 s采集到的圖像中相同位置點的灰度值用前10 s圖像序列對應位置點的像素值代替，也就是保留灌注區域的基礎像素值，消除造影微泡對像素灰度值的影響。提取出的TICs如圖 5所示。前100幀圖像的兩組TICs是重合的，因為此時造影劑還未進入肌瘤，從100幀起造影劑開始進入肌瘤內部，對灌注區域處理后的圖像組的TICs信號強度并未發生劇烈變化，僅有輕微的增強，表明被替換的點是有造影微泡存在的點，即算法準確地識別出了血管灌注區域。
圖5
圖像序列
紅色曲線為運動校正后圖像序列的ROI的TICs，藍色曲線為對灌注區域處理后的圖像序列ROI的TICs
Figure5. Image sequencesred curve： the TICs of ROI that from after correction； blue curve： TICs of ROI from image sequences after processing infusion area
為了進一步驗證提取算法的實用性，由兩位經驗豐富的超聲科醫生觀看造影錄像，再對比由本算法提取出來的血管灌注區域，提取結果和實際灌注區域吻合度高，達到了臨床對灌注區域的識別精度，表明該算法能夠精確地提取出血管灌注區域，具有很好的實用性和準確性。
3 結論
針對子宮肌瘤超聲造影圖像的血管灌注區域的提取，本文提出了一種基于運動校正后的圖像序列閾值提取法。預處理圖像后，首先利用Brox光流法估計出運動場，再根據運動場對相鄰兩幀圖像進行warp校正，接著在腫瘤灌注區域采樣求得分割閾值t,最后根據超聲圖像中富血供與乏血供區域信號強度的變化由閾值t提取出血管灌注區域。結果表明該算法能很好地提取出腫瘤中的血管灌注區域，為臨床對造影區域的識別提供可靠的參考，為進一步超聲造影圖像的定量分析打下基礎。
