本文提出一種基于圖像灰度的關鍵幀提取的門控方法,用于減少血管內超聲(IVUS)圖像序列縱切方向上的運動偽影。首先將IVUS圖像換到極坐標下,通過頻譜分析和濾波技術,提取反映心臟運動的一維信號簇,檢索濾波后信號的極值,尋找關鍵幀,組成門控序列。實驗結果分析表明本文算法快,平均每幀處理時間為17 ms。從IVUS圖像序列的縱向可視圖上觀察,門控序列和原始序列趨勢一致,減少了鋸齒形狀,具有良好的連續性。采集12組臨床IVUS序列[圖像(876±65)幀,血管長度(14.61±1.08) mm],計算門控前后序列的血管容積、管腔容積和平均斑塊負荷。統計實驗結果發現,門控序列血管容積、管腔容積顯著小于原始序列,平均斑塊負荷差異性不顯著,符合臨床診斷需要。血管面積方差和管腔面積方差顯著小于原始序列,表明門控序列較原始序列穩定。
引用本文： 毛海群, 楊豐, 黃錚, 崔凱, 王欣昕. 基于IVUS圖像序列的關鍵幀提取和分析在臨床上的應用. 生物醫學工程學雜志, 2015, 32(4): 892-899. doi: 10.7507/1001-5515.20150159 復制
引言
血管內超聲(intravascular ultrasound，IVUS)是一種有創斷層成像技術。導絲從橈動脈或股動脈穿刺到冠脈血管內，驅動馬達勻速回撤超聲探頭，實時截取血管截面圖像。圖像分辨率高，能準確獲取管壁厚度、管腔形態及動脈粥樣硬化斑塊成分等信息。定量分析IVUS序列血管及管腔面積、容積、最大最小直徑，管腔面積狹窄率，斑塊負荷等診斷指標，為評價冠脈造影難以診斷及隱匿的病變、指導介入治療方案、評價治療效果，序列測定斑塊進展與逆轉等臨床應用提供依據[1]。
由于心臟周期性搏動，引起IVUS序列在短軸平面振蕩，血管尺寸變化，以及冠脈血管部分位置重采樣或未采樣，導致IVUS圖像序列產生運動偽影[2]。偽影在圖像序列中表現非常強烈，在縱切圖中，IVUS圖像序列呈鋸齒狀分布，需要采用計算機輔助的方法，對IVUS圖像序列中存在的運動偽影進行抑制或補償，門控方法是消除運動偽影的方法之一。
門控方法分為：心電圖(electrocardiogram,ECG)門控方法和基于圖像門控方法。ECG門控法又分為在線(on-line)門控[3]和離線(off-line)門控[4]。在線ECG門控法需結合專門的硬件設備，在每個心動周期的R波處采集圖像，導管介入時間長。離線ECG門控同時采集ECG信號和IVUS圖像，醫生對照ECG信號從IVUS序列中標記出門控幀，介入時間減少，但最佳采樣點選擇較難確定，主觀性強。基于圖像的門控法就是完全利用圖像處理技術從IVUS序列中提取類似于ECG的信號，仿照ECG門控方法[5-10]，在每個心動周期的相同階段選擇一幀圖像，組成門控序列。基于圖像的門控法最重要的步驟是：提取IVUS圖像中顯著性特征量作為參考點，參考點隨著IVUS圖像序列變化而形成一個時間變化的序列信號，反映心臟運動的規律。比如用感興趣區域內的平均灰度和絕對灰度差[5]、管腔輪廓的變化[6]作為參考點，使用這類方法的前提需要精準的邊緣檢測方法，保證每幀所提取的顯著性特征量能用于心動周期分析。但目前還沒有精確的IVUS圖像分割方法，限制了該方法的使用。為了降低對特征提取的要求，研究者利用同序列中不同幀圖像之間的非相似性程度來分析心臟運動規律[8-10]，但算法復雜，計算量大。最新的文獻把基于圖像序列的門控法擴展為基于圖像的心動相位(cardiac phase)檢索[11]，其三個步驟為：①選每幀圖像像素點局部均值為參考點，沿序列形成一維信號簇，幅值濾波和光流濾波保留與心臟運動信息相關性很強的一維信號簇；②帶通濾波一維信號簇，組成一個心動周期剖面圖；③檢索心臟相位，通過相位平移找到任意心臟相位對應的門控序列，選擇穩定幀圖像。
所謂的關鍵幀定義為心臟舒張末期對應的IVUS圖像。在心動周期中，舒張末期心室進入等容收縮期，容積不變，心臟運動較慢，超聲探頭和管腔的相對位移小，此時采集到的圖像具有更好的穩定性和一致性[12]。且關鍵幀組成的序列滿足O’Malley提出的門控序列的三個條件[9]：①每個心動周期選取一幀圖像；②門控幀需在心動周期中心臟運動較緩慢的位置處采集；③所有門控幀大致處于同一心臟相位。
本文在文獻[11]的基礎上提出一種改進的基于圖像灰度的關鍵幀門控序列的提取方法。將直角坐標系下的IVUS圖像轉換到極坐標上，去除非成像區域，獲取基于像素灰度變化的一維時間信號。利用頻譜分析技術和濾波技術篩選出反映心臟運動的關鍵點，對關鍵點的灰度變化信號簇進行帶通濾波，只保留反映心臟運動規律信息的信號，檢索帶通濾波信號極值，提取關鍵幀位置。由關鍵幀組成的序列記為關鍵幀門控序列。
1 IVUS關鍵幀提取算法
在IVUS成像采集中，每一幀圖像的像素灰度代表血管的組織結構。冠脈附著在心外膜脂肪內，因此，同一位置像素的灰度變化值代表該點血管組織結構的變化，隱藏了心臟運動變化規律。
1.1 基于IVUS圖像的心動相位檢索原理
文獻[11]指出，IVUS圖像中的像素點可以分為三類：關鍵點、干擾點和背景點。在一組IVUS序列中，取三個像素點(A，B，C)代表三類像素點，它們的灰度均值沿縱切方向變化曲線(longitudinal mean gray,LMG)及LMG信號的幅度譜如圖 1所示。圖 1(a)中有三個像素點：A點(紅色)位于管腔內膜邊緣，B點(藍色)位于管腔接近血管內表面處，C點(品紅色)位于血管附屬組織。圖 1(b)對應像素點的LMG信號，橫坐標為圖像序列幀號，縱坐標為灰度值。圖 1(c)為對應的LMG信號幅度譜，橫坐標為頻率(beats/min)，縱坐標為幅值，右上方框圖為頻率30~140 beats/min范圍內的幅度譜特寫圖。A點的LMG信號呈現出很好的周期性，在頻率50~140 beats/min范圍內存在一個主要頻率成分，是關鍵點，能反映心臟運動規律；B點的LMG信號的幅度譜在頻率50~140 beats/min處分布凌亂，是干擾點；C點的LMG信號分布沒有規律，無有用的頻譜信息，是背景點。
圖1
三個像素點的LMG信號及其LMG信號幅度譜
(a)三個采樣點;(b)采樣點的LMG信號;(c)LMG信號頻域幅度譜
Figure1. LMG signals of three pixels and Fourier amplitude spectrum(a) three sampling pixels; (b) LMG signal of these three pixels; (c) amplitude spectrum of LMG signal
只有關鍵點(如A點)的LMG信號才能反映心臟運動，因此要濾除背景點和干擾點。第一步是背景點的濾除。采用幅值濾波器濾除背景點，統計像素LMG頻域幅值總和，濾除和值較小的點即背景點。第二步是干擾點的濾除。干擾點LMG信號幅度譜的主峰優勢不明顯，而光流濾波器只會選擇主峰頻率成分較大的信號，其表達式為：
| $OF\left( {{f_c}} \right) = \frac{{\hat L\left( {{f_c}} \right)}}{S} - \frac{1}{{\left( {N - 1} \right) \times S}}\sum\limits_{x \in {I_{{f_c}}}/{f_c}} {\hat L\left( x \right)} $ |
式(1)中，fc為心臟估計頻率，Ifc為以fc為中心，半徑為fc/2鄰域，記作U(fc,fc/2)，N為Ifc中的點數，(x)為頻率為x的LMG信號的頻域幅值，，設置閾值OFthr，選擇OF>OFthr的點作為關鍵點。
關鍵點的LMG信號組成一維信號簇，能反映心臟運動。對一維信號簇進行帶通濾波，尋找濾波信號極值，通過相位平移檢索心動相位。
1.2 關鍵幀序列提取方法和步驟
在1.1的基礎上提取關鍵幀門控序列，步驟如下：①將直角坐標下的IVUS圖像轉換到極坐標上，去除非成像區域信息，減少計算量；②分析像素點LMG信號的頻譜特性，設計一種反映幅值、波峰分布的濾波器，篩選關鍵點，提取反映心臟運動規律的一維信號簇；③在頻域中對一維信號簇進行帶通濾波，帶通濾波器的中心頻率為心臟估計頻率，只保留反映心臟運動信息的相關信號；④檢索濾波后的信號極值，獲取心臟舒張末期對應的IVUS圖像即關鍵幀，組成關鍵幀門控序列。
1.2.1 IVUS坐標轉換及像素點采樣
圖 2左邊為一幅IVUS圖像，白色圓域內為成像區域，提供有用信息。以直角坐標系下的圖像中心為圓心(導管中心)，水平向右為0°方向，順時針旋轉生成極坐標圖像。極坐標圖像的橫軸為角度，縱軸為半徑。対一幅512×512的圖像而言，轉換后的計算量減少到原來的35%。
圖2
直角坐標和極坐標下的IVUS圖像
Figure2.
IVUS image in Cartesian and polar coordinates
對極坐標下的圖像進行均值濾波，用像素點的局部灰度均值代替該點的灰度值。對濾波后的圖像再進行非均勻采樣，半徑小的區域采樣點稀疏，半徑大的區域采樣點稠密，保證采樣點在直角坐標系下的分布均勻。同時，極坐標下半徑較小的區域對應導管，半徑較大的區域對應冠脈血管壁外膜以外的附屬結構，這些區域的像素點都是不能反映心臟運動的背景點，為了減少計算量不進行采樣。采樣點在極坐標、直角坐標系下的分布如圖 3(a)上、下圖所示。
圖3
篩選關鍵點
Figure3.
Screen key points
1.2.2 提取反映心臟運動的一維信號簇
由于IVUS圖像中像素點可分為關鍵點、背景點和干擾點。針對三類點LMG信號幅度譜的特點，如圖 1所示，即背景點的LMG信號幅度譜函數在心臟頻率附近處的值很小，干擾點頻譜在心臟頻率附近分布凌亂，主峰優勢不明顯，據此去除這兩類點，篩選關鍵點。
設計一個選擇關鍵點的濾波器。在頻域上，對以心臟頻率為中心的鄰域幅值求和：
| $S\left( {k,{f_c} = } \right)\sum\limits_{x \in {I_{{f_c}}}} {\hat L\left( {k,x} \right)} $ |
光流濾波器只會選擇主峰優勢明顯的信號，表達式為：
| $\begin{array}{l} OF\left( {k,{f_c}} \right) = \frac{{\hat L\left( {k,{f_c}} \right)}}{{S\left( {k,{f_c}} \right)}} - \frac{1}{{\left( {N - 1} \right) \times S\left( {k,{f_c}} \right)}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{x \in {I_{{f_c}}}/{f_c}} {\hat L\left( {k,x} \right)} \end{array}$ |
式(2)、(3)中，fc為心臟估計頻率(在2.2.5節將給出具體討論fc的估計方法)，k為采樣點，(k,x)為點k在頻率x處的LMG信號的頻域幅值，其它參數含義同式(1)。
綜合背景點和干擾點LMG信號的幅度譜特性，設計一種反映幅值、波峰分布的濾波器，選擇關鍵點，表達式為：
| $F\left( {k,{f_c}} \right) = 0.5\frac{{S\left( {k,{f_c}} \right) - {S_{\min }}}}{{{S_{\max }} - {S_{\min }}}} + 0.5\frac{{OF\left( {k,{f_c}} \right) - O{F_{\min }}}}{{O{F_{\max }} - O{F_{\min }}}}$ |
公式(4)中，Smax、Smin、OFmax、OFmin分別為S、OF的最大值、最小值。根據經驗設置閾值Fthr，選擇F>Fthr的點，即為關鍵點，本文選擇F值較大的30個點為關鍵點。
篩選過程如圖 3所示。圖 3(a)為采樣點在極坐標和直角坐標系下的分布圖；按照式(4)計算采樣點的F值，保留F>Fthr的點，即關鍵點。保留下來的點分布如圖 3(b)上部所示，在直角坐標系上的分布如圖 3(b)底部所示。這些點多分布在管壁、管腔和斑塊邊緣，即便部分點在這一幀圖像中沒有分布在邊緣上，在整個IVUS序列上其它的幀圖像上也會體現這一特性。關鍵點的LMG信號分布如圖 4(a)所示，所有關鍵點的LMG信號組成了反映心臟運動的一維信號簇。圖 4(a)表明關鍵點的LMG信號還受到其它因素的干擾，擾亂了其周期性。因此，必須進行帶通濾波，去除這些因素的影響，只保留與心臟信息相關的信號。
圖4
關鍵點的LMG信號及濾波結果圖
(a)關鍵點的LMG信號；(b)濾波信號
Figure4. LMGs and filtered signal of key points(a) LMG signal of a key point; (b) the filtered signal
1.2.3 帶通濾波選擇
臨床上采集到的IVUS序列主要來自心臟疾病的患者，心臟頻率不可能是恒定不變的，同時圖像噪聲、人體呼吸、血管形態結構的變化等因素會干擾一維信號簇的周期性。因此，在檢索關鍵幀之前，必須對關鍵點的LMG信號進行帶通濾波。本文采用以心臟頻率為中心的巴特沃斯帶通濾波器，對關鍵點的LMG信號進行帶通濾波。巴特沃斯帶通濾波器頻域表達式為：
| $B\left( f \right) = {B_{n,b}}\left( f \right) = \frac{1}{{\sqrt {1 + {{\left( {\frac{{{f^2} - f_c^2}}{{bf{f_c}}}} \right)}^{2n}}} }}$ |
式(5)中，fc為心臟估計頻率，n、b為濾波器參數，n與濾波器的衰減速率正相關，b與通帶寬度成正比。
Butterworth帶通濾波器的通帶寬度及衰減快慢對實驗結果存在一定影響。n過小，衰減太慢，濾波后的信號會出現許多的毛刺；n過大，衰減太快，會出現振蕩效應。b越小，濾波后的周期性越好，但由于心率的微小變化，其結果并不十分理想；b越大，濾波后信號一致性不強。通過實驗分析，本文選擇b=0.2，n=3，巴特沃斯濾波器濾波效果如圖 4(b)所示。
1.2.4 關鍵幀門控序列的提取
關鍵點LMG信號的極值對應心臟舒張末期或收縮末期采集到的圖像，通過尋找濾波后的LMG信號極值位置來提取IVUS序列的關鍵幀。極值點對應信號一階導數的過零點。因此，先對濾波后的LMG信號進行三次樣條插值，以保證插值后的新序列的一階導數是平滑的。然后用五點數值微分法求導[13]，導數值趨于零的點就是極值點，極值點所對應的幀為關鍵幀，由此組成關鍵幀門控序列。
1.2.5 心臟頻率估計
篩選關鍵點和設計帶通濾波器需要估計心率值(fc)。心臟運動是產生運動偽影的主要因素，像素點LMG信號的幅度譜在非零頻率附近的最大值所對應的頻率(主峰頻率)作為該點的心臟運動頻率估計值[11]。
通過臨床觀察，在采集IVUS的過程中，患者的心率(beats/min)滿足
| $50 < {f_c} < 140$ |
單純用某一個像素點的主峰頻率估計心率是有偏的，本文采用以下方法估計心率。①計算出所有采樣點的主峰頻率，記為fk，k為采樣點；②取所有滿足公式(6)的fk的均值fh；③重新估計心率，對主峰頻率滿足0.8fh<fk<1.2fh的fk求均值，得到心率估計值fc，用于篩選關鍵點的濾波器[式(4)]設計；④對關鍵點心率重復步驟②，重新得到fc，用于Butterworth帶通濾波器的設計。
2 實驗結果及分析
實驗數據來源于南方醫院7名疑似冠心病患者，取12組IVUS圖像序列。7名疑似冠心病患者的基本信息(性別、年齡)及12組IVUS序列的長度及狹窄情況如表 1所示。序列長度為(876±65)幀，對應的血管長度為(14.61±1.08)mm。數據采集系統是波士頓科學公司的iLabTM Ultrasound Imaging System，超聲探頭為機械式旋轉型，頻率40 MHz，回撤速度0.5 mm/s，采樣頻率為30幀/s。采集到的IVUS圖像像素為512×512，灰度為256階灰度。實驗平臺采用個人計算機，Matlab2011a編程。
表1
7名疑似冠心病患者的基本信息及12組IVUS序列的長度及狹窄情況
Table1.
Basic information of 7 suspicious patients with coronary heart disease and 12 IVUS sets
2.1 關鍵幀提取
采用本文方法提取IVUS圖像序列的關鍵幀，算法較快，表 1最后一列列出了每組序列的處理時間，平均每幀圖像的處理時間為(17.46±1.07) ms。圖 5為一組IVUS圖像序列(1 000幀)的縱切圖，管壁呈鋸齒狀分布，圖中白色豎直線代表關鍵幀。圖 5表明，關鍵幀主要分布在鋸齒狀波形的峰值處。
圖5
一組IVUS縱切圖及關鍵幀
Figure5.
A non-gated IVUS sequence longitudinal cut and its key frames
圖 6為一組IVUS序列(1 000幀)在不同的縱切方向上(0°、45°、90°、145°)的縱切圖。圖 7為4組IVUS圖像序列的關鍵幀檢索結果，序列長度均為1 000幀。 視覺上看，關鍵幀門控序列管壁、斑塊與原始序列具有相似的趨勢，血管內腔和管壁輪廓和斑塊比原始序列要平滑，連續性更好，抑制了鋸齒狀振蕩，具有更好的視覺效果。
圖6
同一組序列不同方向(0°，45°，90°，145°)的原始和關鍵幀門控序列縱視圖
Figure6.
An IVUS set longitudinal cuts at different angles (0°，45°，90°，145°) and their key-frames gated cuts
圖7
4組IVUS序列原始序列及其門控序列
Figure7.
Four original IVUS sequences longitudinal cuts and their key-frames gated sequences
2.2 參數分析
2.2.1 IVUS圖像序列參數計算方法
在IVUS冠心病診療中，常用到的診斷指標有：血管壁，管腔面積和容積，斑塊面積，斑塊容積，管腔最大最小直徑，管腔面積狹窄率，斑塊負荷等。血管容積、管腔容積、斑塊負荷等是最重要的參數，這些參數的定量計算需要檢測管腔/內膜(lumen-intima)界面和外彈力膜(external elastic membrane,EEM)界面這兩個主要分界面，其中，EEM界面近似為血管壁，lumen-intima界面近似為管腔邊界。
EEM界面近似為橢圓[14]，在2-D IVUS圖像中標記EEM的八個邊界點，用最小二乘法擬合橢圓的一般方程[15]。下面將具體介紹最小二乘法擬合橢圓方程的方法，橢圓的一般方程為：
| $a{x^2} + bxy + c{y^2} + dx + ey + f = 0$ |
令α=[a,b,c,d,e,f]，X=[x2,xy,y2,x,y,1]，故橢圓一般方程可以表示為:
| $f\left( {\boldsymbol{\alpha} ,\boldsymbol{X}} \right) = \boldsymbol{\alpha} {\boldsymbol{X}^T} = a{x^2} + bxy + c{y^2} + dx + ey + f = 0$ |
假設(xi,yi)是標記出的EEM邊界點，根據最小二乘原理，曲線擬合問題可以通過代數距離平方和最小化實現，代價函數為
| $Cost\left( \boldsymbol{\alpha} \right) = \sum\limits_{i = 1}^8 {{{\left( {ax_i^2 + b{x_i}{y_i} + cy_i^2 + d{x_i} + e{y_i} + f} \right)}^2}} $ |
由極值原理，要使得Cost(α)最小，必有：
| $\frac{{\partial Cost}}{{\partial a}} = \frac{{\partial Cost}}{{\partial b}} = \frac{{\partial Cost}}{{\partial c}} = \frac{{\partial Cost}}{{\partial d}} = \frac{{\partial Cost}}{{\partial e}} = \frac{{\partial Cost}}{{\partial f}} = 0$ |
由此可得到橢圓方程的系數，從而得到橢圓的一般方程，根據橢圓一般方程計算橢圓長軸、短軸、中心及面積，據此確定血管面積和最大最小直徑。
目前還沒有確定的模型模擬lumen-intima界面，手動標記16個邊界點，樣條插值擬合管腔邊界。管腔中心定義為: ，[μi,vi]為擬合到的邊界點坐標，N為管腔/內膜邊界點數[16]。確定管腔面積、最大最小直徑。結合兩界面參數得到斑塊負荷等血管參數。
原始序列和門控序列血管容積分別定義為：
| ${V_{oV}} = \sum\limits_{i = 1}^n {0.5} \times \left[ {{S_{oV}}\left( i \right) + {S_{oV}}\left( {i + 1} \right)} \right] \times d$ |
| ${V_{gV}} = \sum\limits_{i = 1}^n {0.5} \times \left[ {{S_{gV}}\left( i \right) + {S_{gV}}\left( {i + 1} \right)} \right] \times h\left( i \right) \times d$ |
式(11)、(12)中，N、n分別是原始序列及門控序列幀數，SoV(i)、SgV(i)分別是原始序列和門控序列第i幀圖像的血管橫截面積，d為一幀圖像對應的血管厚度，d=1/60 mm，h(i)為相鄰門控幀之間的真實幀間距。管腔容積的計算類似。
斑塊負荷定義為為：
| ${B_p} = \left( {{S_\upsilon } - {S_l}} \right)/{S_\upsilon }$ |
式(13)中，Sv、Sl分別為血管橫截面積，管腔橫截面積。
圖 8為一組含有斑塊的IVUS圖像序列(1 000幀)原始序列血管面積及管腔面積，門控序列血管面積及管腔面積的變化曲線。門控序列面積變化比原始序列平緩，保留了原始序列中的信息，具有更加直觀的視覺效果。
圖8
原始序列和關鍵幀門控序列血管及管腔面積變化曲線
Figure8.
Vessel and lumen cross-sectional area of an original and key-frames gated set
2.2.2 原始序列及門控序列參數比較
在臨床醫生的指導下，人工畫出12組IVUS圖像序列的EEM界面和lumen-intima界面。根據2.2.1 節給出方法，計算門控前后序列的血管容積、管腔容積、平均斑塊負荷，結果如表 2所示，表 2中IVUS序列編號與表 1序列編號是一一對應的。
表2
12組IVUS圖像原始序列及關鍵幀門控序列的血管參數和分析統計
Table2.
Vascular parameters and statistic analysis on 12 IVUS non-gated vs. key-frames gated sequences
正常情況下，冠脈血管的管徑在心臟收縮末期達到最大值，心臟舒張末期冠脈管徑相對較小[1]。關鍵幀對應心臟舒張末期的圖像，其血管和管腔橫截面積相對較小，圖 8也說明了這一點。因此，門控序列的血管容積、管腔容積應小于原始序列。配對資料的符號秩和檢驗適用于配對設計的資料，使用對資料的符號秩和檢驗檢驗原始序列及門控序列血管容積、管腔容積、平均斑塊負荷這三個參數的差異性是否顯著。檢驗結果如表 2最后一行所示，配對秩和檢驗的置信水平為α=0.05。檢驗結果顯示，門控序列的管腔容積、血管容積顯著小于原始序列，符合實際；門控前后的平均斑塊負荷差異性不顯著，表明門控序列保留了原始序列中管腔面積與血管面積之間的相關信息，符合臨床診斷需要。
IVUS原始序列和門控序列的血管橫截面積、管腔橫截面積的k檢驗結果表明它們都符合正態分布。對正態分布下每一組原始序列和門控序列的血管橫截面積方差進行F檢驗。k檢驗和F檢驗的置信水平均為α=0.05。檢驗結果表明門控序列的血管面積方差顯著小于原始序列，同樣，門控序列的管腔面積方差顯著小于原始序列，表明門控序列血管橫截面積和管腔橫截面積都比原始序列穩定。
3 結論
本文提出了一種簡單、基本實時的基于IVUS圖像序列灰度的門控方法。該方法基于圖像像素灰度值和組織結構特性的一致性，用空間像素點的LMG信號作為圖像中顯著性特征量，提取反映心臟運動規律的一維信號。定量分析了門控前后序列的各個參數之間的關系，門控序列包含了原始序列中的關鍵信息，減少了臨床血管參數的分析量，具有和原始序列等同甚至更好的臨床參考價值，相比于原始序列，門控序列具有更好的視覺效果。該方法也存在一定的局限性：冠脈血管結構多變，分支、毗嶺結構、斑塊、血栓等結構可能交替出現，對于這類序列，本文算法難以準確檢索到其關鍵幀門控序列。其次，在設計帶通濾波器的時候中心頻率是恒定的，對于心率不齊的患者，難以準確提取其門控序列。
引言
血管內超聲(intravascular ultrasound，IVUS)是一種有創斷層成像技術。導絲從橈動脈或股動脈穿刺到冠脈血管內，驅動馬達勻速回撤超聲探頭，實時截取血管截面圖像。圖像分辨率高，能準確獲取管壁厚度、管腔形態及動脈粥樣硬化斑塊成分等信息。定量分析IVUS序列血管及管腔面積、容積、最大最小直徑，管腔面積狹窄率，斑塊負荷等診斷指標，為評價冠脈造影難以診斷及隱匿的病變、指導介入治療方案、評價治療效果，序列測定斑塊進展與逆轉等臨床應用提供依據[1]。
由于心臟周期性搏動，引起IVUS序列在短軸平面振蕩，血管尺寸變化，以及冠脈血管部分位置重采樣或未采樣，導致IVUS圖像序列產生運動偽影[2]。偽影在圖像序列中表現非常強烈，在縱切圖中，IVUS圖像序列呈鋸齒狀分布，需要采用計算機輔助的方法，對IVUS圖像序列中存在的運動偽影進行抑制或補償，門控方法是消除運動偽影的方法之一。
門控方法分為：心電圖(electrocardiogram,ECG)門控方法和基于圖像門控方法。ECG門控法又分為在線(on-line)門控[3]和離線(off-line)門控[4]。在線ECG門控法需結合專門的硬件設備，在每個心動周期的R波處采集圖像，導管介入時間長。離線ECG門控同時采集ECG信號和IVUS圖像，醫生對照ECG信號從IVUS序列中標記出門控幀，介入時間減少，但最佳采樣點選擇較難確定，主觀性強。基于圖像的門控法就是完全利用圖像處理技術從IVUS序列中提取類似于ECG的信號，仿照ECG門控方法[5-10]，在每個心動周期的相同階段選擇一幀圖像，組成門控序列。基于圖像的門控法最重要的步驟是：提取IVUS圖像中顯著性特征量作為參考點，參考點隨著IVUS圖像序列變化而形成一個時間變化的序列信號，反映心臟運動的規律。比如用感興趣區域內的平均灰度和絕對灰度差[5]、管腔輪廓的變化[6]作為參考點，使用這類方法的前提需要精準的邊緣檢測方法，保證每幀所提取的顯著性特征量能用于心動周期分析。但目前還沒有精確的IVUS圖像分割方法，限制了該方法的使用。為了降低對特征提取的要求，研究者利用同序列中不同幀圖像之間的非相似性程度來分析心臟運動規律[8-10]，但算法復雜，計算量大。最新的文獻把基于圖像序列的門控法擴展為基于圖像的心動相位(cardiac phase)檢索[11]，其三個步驟為：①選每幀圖像像素點局部均值為參考點，沿序列形成一維信號簇，幅值濾波和光流濾波保留與心臟運動信息相關性很強的一維信號簇；②帶通濾波一維信號簇，組成一個心動周期剖面圖；③檢索心臟相位，通過相位平移找到任意心臟相位對應的門控序列，選擇穩定幀圖像。
所謂的關鍵幀定義為心臟舒張末期對應的IVUS圖像。在心動周期中，舒張末期心室進入等容收縮期，容積不變，心臟運動較慢，超聲探頭和管腔的相對位移小，此時采集到的圖像具有更好的穩定性和一致性[12]。且關鍵幀組成的序列滿足O’Malley提出的門控序列的三個條件[9]：①每個心動周期選取一幀圖像；②門控幀需在心動周期中心臟運動較緩慢的位置處采集；③所有門控幀大致處于同一心臟相位。
本文在文獻[11]的基礎上提出一種改進的基于圖像灰度的關鍵幀門控序列的提取方法。將直角坐標系下的IVUS圖像轉換到極坐標上，去除非成像區域，獲取基于像素灰度變化的一維時間信號。利用頻譜分析技術和濾波技術篩選出反映心臟運動的關鍵點，對關鍵點的灰度變化信號簇進行帶通濾波，只保留反映心臟運動規律信息的信號，檢索帶通濾波信號極值，提取關鍵幀位置。由關鍵幀組成的序列記為關鍵幀門控序列。
1 IVUS關鍵幀提取算法
在IVUS成像采集中，每一幀圖像的像素灰度代表血管的組織結構。冠脈附著在心外膜脂肪內，因此，同一位置像素的灰度變化值代表該點血管組織結構的變化，隱藏了心臟運動變化規律。
1.1 基于IVUS圖像的心動相位檢索原理
文獻[11]指出，IVUS圖像中的像素點可以分為三類：關鍵點、干擾點和背景點。在一組IVUS序列中，取三個像素點(A，B，C)代表三類像素點，它們的灰度均值沿縱切方向變化曲線(longitudinal mean gray,LMG)及LMG信號的幅度譜如圖 1所示。圖 1(a)中有三個像素點：A點(紅色)位于管腔內膜邊緣，B點(藍色)位于管腔接近血管內表面處，C點(品紅色)位于血管附屬組織。圖 1(b)對應像素點的LMG信號，橫坐標為圖像序列幀號，縱坐標為灰度值。圖 1(c)為對應的LMG信號幅度譜，橫坐標為頻率(beats/min)，縱坐標為幅值，右上方框圖為頻率30~140 beats/min范圍內的幅度譜特寫圖。A點的LMG信號呈現出很好的周期性，在頻率50~140 beats/min范圍內存在一個主要頻率成分，是關鍵點，能反映心臟運動規律；B點的LMG信號的幅度譜在頻率50~140 beats/min處分布凌亂，是干擾點；C點的LMG信號分布沒有規律，無有用的頻譜信息，是背景點。
圖1
三個像素點的LMG信號及其LMG信號幅度譜
(a)三個采樣點;(b)采樣點的LMG信號;(c)LMG信號頻域幅度譜
Figure1. LMG signals of three pixels and Fourier amplitude spectrum(a) three sampling pixels; (b) LMG signal of these three pixels; (c) amplitude spectrum of LMG signal
只有關鍵點(如A點)的LMG信號才能反映心臟運動，因此要濾除背景點和干擾點。第一步是背景點的濾除。采用幅值濾波器濾除背景點，統計像素LMG頻域幅值總和，濾除和值較小的點即背景點。第二步是干擾點的濾除。干擾點LMG信號幅度譜的主峰優勢不明顯，而光流濾波器只會選擇主峰頻率成分較大的信號，其表達式為：
| $OF\left( {{f_c}} \right) = \frac{{\hat L\left( {{f_c}} \right)}}{S} - \frac{1}{{\left( {N - 1} \right) \times S}}\sum\limits_{x \in {I_{{f_c}}}/{f_c}} {\hat L\left( x \right)} $ |
式(1)中，fc為心臟估計頻率，Ifc為以fc為中心，半徑為fc/2鄰域，記作U(fc,fc/2)，N為Ifc中的點數，(x)為頻率為x的LMG信號的頻域幅值，，設置閾值OFthr，選擇OF>OFthr的點作為關鍵點。
關鍵點的LMG信號組成一維信號簇，能反映心臟運動。對一維信號簇進行帶通濾波，尋找濾波信號極值，通過相位平移檢索心動相位。
1.2 關鍵幀序列提取方法和步驟
在1.1的基礎上提取關鍵幀門控序列，步驟如下：①將直角坐標下的IVUS圖像轉換到極坐標上，去除非成像區域信息，減少計算量；②分析像素點LMG信號的頻譜特性，設計一種反映幅值、波峰分布的濾波器，篩選關鍵點，提取反映心臟運動規律的一維信號簇；③在頻域中對一維信號簇進行帶通濾波，帶通濾波器的中心頻率為心臟估計頻率，只保留反映心臟運動信息的相關信號；④檢索濾波后的信號極值，獲取心臟舒張末期對應的IVUS圖像即關鍵幀，組成關鍵幀門控序列。
1.2.1 IVUS坐標轉換及像素點采樣
圖 2左邊為一幅IVUS圖像，白色圓域內為成像區域，提供有用信息。以直角坐標系下的圖像中心為圓心(導管中心)，水平向右為0°方向，順時針旋轉生成極坐標圖像。極坐標圖像的橫軸為角度，縱軸為半徑。対一幅512×512的圖像而言，轉換后的計算量減少到原來的35%。
圖2
直角坐標和極坐標下的IVUS圖像
Figure2.
IVUS image in Cartesian and polar coordinates
對極坐標下的圖像進行均值濾波，用像素點的局部灰度均值代替該點的灰度值。對濾波后的圖像再進行非均勻采樣，半徑小的區域采樣點稀疏，半徑大的區域采樣點稠密，保證采樣點在直角坐標系下的分布均勻。同時，極坐標下半徑較小的區域對應導管，半徑較大的區域對應冠脈血管壁外膜以外的附屬結構，這些區域的像素點都是不能反映心臟運動的背景點，為了減少計算量不進行采樣。采樣點在極坐標、直角坐標系下的分布如圖 3(a)上、下圖所示。
圖3
篩選關鍵點
Figure3.
Screen key points
1.2.2 提取反映心臟運動的一維信號簇
由于IVUS圖像中像素點可分為關鍵點、背景點和干擾點。針對三類點LMG信號幅度譜的特點，如圖 1所示，即背景點的LMG信號幅度譜函數在心臟頻率附近處的值很小，干擾點頻譜在心臟頻率附近分布凌亂，主峰優勢不明顯，據此去除這兩類點，篩選關鍵點。
設計一個選擇關鍵點的濾波器。在頻域上，對以心臟頻率為中心的鄰域幅值求和：
| $S\left( {k,{f_c} = } \right)\sum\limits_{x \in {I_{{f_c}}}} {\hat L\left( {k,x} \right)} $ |
光流濾波器只會選擇主峰優勢明顯的信號，表達式為：
| $\begin{array}{l} OF\left( {k,{f_c}} \right) = \frac{{\hat L\left( {k,{f_c}} \right)}}{{S\left( {k,{f_c}} \right)}} - \frac{1}{{\left( {N - 1} \right) \times S\left( {k,{f_c}} \right)}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{x \in {I_{{f_c}}}/{f_c}} {\hat L\left( {k,x} \right)} \end{array}$ |
式(2)、(3)中，fc為心臟估計頻率(在2.2.5節將給出具體討論fc的估計方法)，k為采樣點，(k,x)為點k在頻率x處的LMG信號的頻域幅值，其它參數含義同式(1)。
綜合背景點和干擾點LMG信號的幅度譜特性，設計一種反映幅值、波峰分布的濾波器，選擇關鍵點，表達式為：
| $F\left( {k,{f_c}} \right) = 0.5\frac{{S\left( {k,{f_c}} \right) - {S_{\min }}}}{{{S_{\max }} - {S_{\min }}}} + 0.5\frac{{OF\left( {k,{f_c}} \right) - O{F_{\min }}}}{{O{F_{\max }} - O{F_{\min }}}}$ |
公式(4)中，Smax、Smin、OFmax、OFmin分別為S、OF的最大值、最小值。根據經驗設置閾值Fthr，選擇F>Fthr的點，即為關鍵點，本文選擇F值較大的30個點為關鍵點。
篩選過程如圖 3所示。圖 3(a)為采樣點在極坐標和直角坐標系下的分布圖；按照式(4)計算采樣點的F值，保留F>Fthr的點，即關鍵點。保留下來的點分布如圖 3(b)上部所示，在直角坐標系上的分布如圖 3(b)底部所示。這些點多分布在管壁、管腔和斑塊邊緣，即便部分點在這一幀圖像中沒有分布在邊緣上，在整個IVUS序列上其它的幀圖像上也會體現這一特性。關鍵點的LMG信號分布如圖 4(a)所示，所有關鍵點的LMG信號組成了反映心臟運動的一維信號簇。圖 4(a)表明關鍵點的LMG信號還受到其它因素的干擾，擾亂了其周期性。因此，必須進行帶通濾波，去除這些因素的影響，只保留與心臟信息相關的信號。
圖4
關鍵點的LMG信號及濾波結果圖
(a)關鍵點的LMG信號；(b)濾波信號
Figure4. LMGs and filtered signal of key points(a) LMG signal of a key point; (b) the filtered signal
1.2.3 帶通濾波選擇
臨床上采集到的IVUS序列主要來自心臟疾病的患者，心臟頻率不可能是恒定不變的，同時圖像噪聲、人體呼吸、血管形態結構的變化等因素會干擾一維信號簇的周期性。因此，在檢索關鍵幀之前，必須對關鍵點的LMG信號進行帶通濾波。本文采用以心臟頻率為中心的巴特沃斯帶通濾波器，對關鍵點的LMG信號進行帶通濾波。巴特沃斯帶通濾波器頻域表達式為：
| $B\left( f \right) = {B_{n,b}}\left( f \right) = \frac{1}{{\sqrt {1 + {{\left( {\frac{{{f^2} - f_c^2}}{{bf{f_c}}}} \right)}^{2n}}} }}$ |
式(5)中，fc為心臟估計頻率，n、b為濾波器參數，n與濾波器的衰減速率正相關，b與通帶寬度成正比。
Butterworth帶通濾波器的通帶寬度及衰減快慢對實驗結果存在一定影響。n過小，衰減太慢，濾波后的信號會出現許多的毛刺；n過大，衰減太快，會出現振蕩效應。b越小，濾波后的周期性越好，但由于心率的微小變化，其結果并不十分理想；b越大，濾波后信號一致性不強。通過實驗分析，本文選擇b=0.2，n=3，巴特沃斯濾波器濾波效果如圖 4(b)所示。
1.2.4 關鍵幀門控序列的提取
關鍵點LMG信號的極值對應心臟舒張末期或收縮末期采集到的圖像，通過尋找濾波后的LMG信號極值位置來提取IVUS序列的關鍵幀。極值點對應信號一階導數的過零點。因此，先對濾波后的LMG信號進行三次樣條插值，以保證插值后的新序列的一階導數是平滑的。然后用五點數值微分法求導[13]，導數值趨于零的點就是極值點，極值點所對應的幀為關鍵幀，由此組成關鍵幀門控序列。
1.2.5 心臟頻率估計
篩選關鍵點和設計帶通濾波器需要估計心率值(fc)。心臟運動是產生運動偽影的主要因素，像素點LMG信號的幅度譜在非零頻率附近的最大值所對應的頻率(主峰頻率)作為該點的心臟運動頻率估計值[11]。
通過臨床觀察，在采集IVUS的過程中，患者的心率(beats/min)滿足
| $50 < {f_c} < 140$ |
單純用某一個像素點的主峰頻率估計心率是有偏的，本文采用以下方法估計心率。①計算出所有采樣點的主峰頻率，記為fk，k為采樣點；②取所有滿足公式(6)的fk的均值fh；③重新估計心率，對主峰頻率滿足0.8fh<fk<1.2fh的fk求均值，得到心率估計值fc，用于篩選關鍵點的濾波器[式(4)]設計；④對關鍵點心率重復步驟②，重新得到fc，用于Butterworth帶通濾波器的設計。
2 實驗結果及分析
實驗數據來源于南方醫院7名疑似冠心病患者，取12組IVUS圖像序列。7名疑似冠心病患者的基本信息(性別、年齡)及12組IVUS序列的長度及狹窄情況如表 1所示。序列長度為(876±65)幀，對應的血管長度為(14.61±1.08)mm。數據采集系統是波士頓科學公司的iLabTM Ultrasound Imaging System，超聲探頭為機械式旋轉型，頻率40 MHz，回撤速度0.5 mm/s，采樣頻率為30幀/s。采集到的IVUS圖像像素為512×512，灰度為256階灰度。實驗平臺采用個人計算機，Matlab2011a編程。
表1
7名疑似冠心病患者的基本信息及12組IVUS序列的長度及狹窄情況
Table1.
Basic information of 7 suspicious patients with coronary heart disease and 12 IVUS sets
2.1 關鍵幀提取
采用本文方法提取IVUS圖像序列的關鍵幀，算法較快，表 1最后一列列出了每組序列的處理時間，平均每幀圖像的處理時間為(17.46±1.07) ms。圖 5為一組IVUS圖像序列(1 000幀)的縱切圖，管壁呈鋸齒狀分布，圖中白色豎直線代表關鍵幀。圖 5表明，關鍵幀主要分布在鋸齒狀波形的峰值處。
圖5
一組IVUS縱切圖及關鍵幀
Figure5.
A non-gated IVUS sequence longitudinal cut and its key frames
圖 6為一組IVUS序列(1 000幀)在不同的縱切方向上(0°、45°、90°、145°)的縱切圖。圖 7為4組IVUS圖像序列的關鍵幀檢索結果，序列長度均為1 000幀。 視覺上看，關鍵幀門控序列管壁、斑塊與原始序列具有相似的趨勢，血管內腔和管壁輪廓和斑塊比原始序列要平滑，連續性更好，抑制了鋸齒狀振蕩，具有更好的視覺效果。
圖6
同一組序列不同方向(0°，45°，90°，145°)的原始和關鍵幀門控序列縱視圖
Figure6.
An IVUS set longitudinal cuts at different angles (0°，45°，90°，145°) and their key-frames gated cuts
圖7
4組IVUS序列原始序列及其門控序列
Figure7.
Four original IVUS sequences longitudinal cuts and their key-frames gated sequences
2.2 參數分析
2.2.1 IVUS圖像序列參數計算方法
在IVUS冠心病診療中，常用到的診斷指標有：血管壁，管腔面積和容積，斑塊面積，斑塊容積，管腔最大最小直徑，管腔面積狹窄率，斑塊負荷等。血管容積、管腔容積、斑塊負荷等是最重要的參數，這些參數的定量計算需要檢測管腔/內膜(lumen-intima)界面和外彈力膜(external elastic membrane,EEM)界面這兩個主要分界面，其中，EEM界面近似為血管壁，lumen-intima界面近似為管腔邊界。
EEM界面近似為橢圓[14]，在2-D IVUS圖像中標記EEM的八個邊界點，用最小二乘法擬合橢圓的一般方程[15]。下面將具體介紹最小二乘法擬合橢圓方程的方法，橢圓的一般方程為：
| $a{x^2} + bxy + c{y^2} + dx + ey + f = 0$ |
令α=[a,b,c,d,e,f]，X=[x2,xy,y2,x,y,1]，故橢圓一般方程可以表示為:
| $f\left( {\boldsymbol{\alpha} ,\boldsymbol{X}} \right) = \boldsymbol{\alpha} {\boldsymbol{X}^T} = a{x^2} + bxy + c{y^2} + dx + ey + f = 0$ |
假設(xi,yi)是標記出的EEM邊界點，根據最小二乘原理，曲線擬合問題可以通過代數距離平方和最小化實現，代價函數為
| $Cost\left( \boldsymbol{\alpha} \right) = \sum\limits_{i = 1}^8 {{{\left( {ax_i^2 + b{x_i}{y_i} + cy_i^2 + d{x_i} + e{y_i} + f} \right)}^2}} $ |
由極值原理，要使得Cost(α)最小，必有：
| $\frac{{\partial Cost}}{{\partial a}} = \frac{{\partial Cost}}{{\partial b}} = \frac{{\partial Cost}}{{\partial c}} = \frac{{\partial Cost}}{{\partial d}} = \frac{{\partial Cost}}{{\partial e}} = \frac{{\partial Cost}}{{\partial f}} = 0$ |
由此可得到橢圓方程的系數，從而得到橢圓的一般方程，根據橢圓一般方程計算橢圓長軸、短軸、中心及面積，據此確定血管面積和最大最小直徑。
目前還沒有確定的模型模擬lumen-intima界面，手動標記16個邊界點，樣條插值擬合管腔邊界。管腔中心定義為: ，[μi,vi]為擬合到的邊界點坐標，N為管腔/內膜邊界點數[16]。確定管腔面積、最大最小直徑。結合兩界面參數得到斑塊負荷等血管參數。
原始序列和門控序列血管容積分別定義為：
| ${V_{oV}} = \sum\limits_{i = 1}^n {0.5} \times \left[ {{S_{oV}}\left( i \right) + {S_{oV}}\left( {i + 1} \right)} \right] \times d$ |
| ${V_{gV}} = \sum\limits_{i = 1}^n {0.5} \times \left[ {{S_{gV}}\left( i \right) + {S_{gV}}\left( {i + 1} \right)} \right] \times h\left( i \right) \times d$ |
式(11)、(12)中，N、n分別是原始序列及門控序列幀數，SoV(i)、SgV(i)分別是原始序列和門控序列第i幀圖像的血管橫截面積，d為一幀圖像對應的血管厚度，d=1/60 mm，h(i)為相鄰門控幀之間的真實幀間距。管腔容積的計算類似。
斑塊負荷定義為為：
| ${B_p} = \left( {{S_\upsilon } - {S_l}} \right)/{S_\upsilon }$ |
式(13)中，Sv、Sl分別為血管橫截面積，管腔橫截面積。
圖 8為一組含有斑塊的IVUS圖像序列(1 000幀)原始序列血管面積及管腔面積，門控序列血管面積及管腔面積的變化曲線。門控序列面積變化比原始序列平緩，保留了原始序列中的信息，具有更加直觀的視覺效果。
圖8
原始序列和關鍵幀門控序列血管及管腔面積變化曲線
Figure8.
Vessel and lumen cross-sectional area of an original and key-frames gated set
2.2.2 原始序列及門控序列參數比較
在臨床醫生的指導下，人工畫出12組IVUS圖像序列的EEM界面和lumen-intima界面。根據2.2.1 節給出方法，計算門控前后序列的血管容積、管腔容積、平均斑塊負荷，結果如表 2所示，表 2中IVUS序列編號與表 1序列編號是一一對應的。
表2
12組IVUS圖像原始序列及關鍵幀門控序列的血管參數和分析統計
Table2.
Vascular parameters and statistic analysis on 12 IVUS non-gated vs. key-frames gated sequences
正常情況下，冠脈血管的管徑在心臟收縮末期達到最大值，心臟舒張末期冠脈管徑相對較小[1]。關鍵幀對應心臟舒張末期的圖像，其血管和管腔橫截面積相對較小，圖 8也說明了這一點。因此，門控序列的血管容積、管腔容積應小于原始序列。配對資料的符號秩和檢驗適用于配對設計的資料，使用對資料的符號秩和檢驗檢驗原始序列及門控序列血管容積、管腔容積、平均斑塊負荷這三個參數的差異性是否顯著。檢驗結果如表 2最后一行所示，配對秩和檢驗的置信水平為α=0.05。檢驗結果顯示，門控序列的管腔容積、血管容積顯著小于原始序列，符合實際；門控前后的平均斑塊負荷差異性不顯著，表明門控序列保留了原始序列中管腔面積與血管面積之間的相關信息，符合臨床診斷需要。
IVUS原始序列和門控序列的血管橫截面積、管腔橫截面積的k檢驗結果表明它們都符合正態分布。對正態分布下每一組原始序列和門控序列的血管橫截面積方差進行F檢驗。k檢驗和F檢驗的置信水平均為α=0.05。檢驗結果表明門控序列的血管面積方差顯著小于原始序列，同樣，門控序列的管腔面積方差顯著小于原始序列，表明門控序列血管橫截面積和管腔橫截面積都比原始序列穩定。
3 結論
本文提出了一種簡單、基本實時的基于IVUS圖像序列灰度的門控方法。該方法基于圖像像素灰度值和組織結構特性的一致性，用空間像素點的LMG信號作為圖像中顯著性特征量，提取反映心臟運動規律的一維信號。定量分析了門控前后序列的各個參數之間的關系，門控序列包含了原始序列中的關鍵信息，減少了臨床血管參數的分析量，具有和原始序列等同甚至更好的臨床參考價值，相比于原始序列，門控序列具有更好的視覺效果。該方法也存在一定的局限性：冠脈血管結構多變，分支、毗嶺結構、斑塊、血栓等結構可能交替出現，對于這類序列，本文算法難以準確檢索到其關鍵幀門控序列。其次，在設計帶通濾波器的時候中心頻率是恒定的，對于心率不齊的患者，難以準確提取其門控序列。
