(北京工業大學 信息學部,北京 100124)
隨著年齡的增長,心血管事件如冠狀動脈疾病或動脈粥樣硬化,常常引起生物組織硬度變化和動脈粥樣硬化斑塊的產生,最終導致中風或心肌梗死。最初,將安裝超聲儀的導管進入血管內獲得IVUS成像,來識別動脈壁形態,纖維斑塊位置和厚度以及內腔面積,改善心血管疾病(如動脈粥樣硬化)的診斷。因此,其臨床應用局限于輔助性的超聲波檢查工具或精湛的血管內手術。
然而,彈性成像是一種非常有前景的方法來表征血管壁的力學性能。分析血管彈性對預防斑塊硬化、易損斑塊的早期發現以及指導動脈治療手術具有十分重要的意義[1]。生物力學研究表明[3],動脈粥樣硬化斑塊某個位置的應力超過局部閾值時易發生破裂。但是,如何準確評價血管斑塊彈性及預測斑塊破裂風險的方法仍然缺乏[4]。
非剛性圖像配準一直是醫學圖像處理領域研究的熱點[5],可以用來計算血管壁變形,研究斑塊破裂風險,從而為預防和治療心血管疾病提供依據[6-7]。LiangYun等人[7]使用有限元網格和三次Bspline基函數作為配準模型,采用列文伯格-馬夸爾特(LM)法方法最小化配準函數。然而,三次Bsplines的離散初始化過分平滑應變值,使測量的應變不準確。Richards等人[9-10]使用半范數正則化配準模型,采用牛頓下降法最小化配準函數。但位移較大時,變形效果不太真實。Tang Z[11]等人在配準能模型引入平衡系數,一定程度緩解了配準精度較低問題,但計算速度較慢。
由于處理IVUS圖像本身不足以解決由于高噪聲造成的位移模糊問題,因此使用正則化方法為整合這些約束信息提供了必要的理論框架。由于血管內外膜運動的復雜性和空間異質性,全局參數正則化模型是不合適的。而光流方法是一個更好的模型,它允許估計局部變形。因此,為了進一度提高位移估計的配準精度和計算速度,本文提出一種結合光流的非剛性配準模型,使用改進的多尺度方法最小化配準函數,得到了較好的位移結果。
非剛性圖像配準被認為是一種優化問題,其目標是最小化配準函數。 在本文中,所提出的配準函數包括表示前后兩幀IVUS圖像的相似性項,以及結合組織位移平滑約束增加位移魯棒性的加權項。
經典的配準模型[12],如式(1):
E=Sim+λReg
(1)
E表示圖像配準函數,Sim表示圖像配準項,用來評價前后兩幀圖像的相關性。Reg表示正則項,用來平滑位移函數。λ表示正則項比重。
其中,Sim采用平方和差(SSD)準則[12],,如式(2):
(2)
I(x,y)代表前一幀IVUS圖像,J(x,y)代表后一幀IVUS圖像。D表示圖像I和圖像J的變換矩陣,即待求解的位移場。Ω表示IVUS圖像中計算位移的區域。
Reg采用灰度平滑約束和光流梯度約束,用來提高位移估計的魯棒性,如式(3):
Reg=(▽I(x,y)-▽J(x,y))2+Q(▽D:▽D)
(3)
▽I(x,y)表示前一幀IVUS圖像的灰度梯度,▽J(x,y)表示后一幀IVUS圖像的灰度梯度,(▽I(x,y)-▽J(x,y))2表示灰度平滑約束。▽D表示光流梯度,:符號表示矩陣內積,(▽D:▽D)表示光流梯度約束。Q表示空間權重函數,用于調節光流梯度約束,本文以血管內膜距血管中心的徑向距離R(x,y)表示即Q=R4。
將式(2)和式(3)代入式(1)中得到最終配準模型,如式(4):

(4)
最小化能量函數是一個經典的優化問題,能量函數是非凸的,可能包含幾個局部最小值。其次,最小化計算時間是一個重要的考慮因素。如模擬退火等隨機算法在理論上向能量函數的全局最小值收斂。然而,隨機松弛算法需要大量的迭代非常緩慢;迭代條件模式(ICM)算法等確定性算法會更快地收斂,但可以能陷入能量函數的局部極小值中。為了避免選擇局部最小值,同時減少迭代時間,多尺度方法通過在更精細的配置空間內搜索提供了這種可能性。
多尺度迭代方法[13-14]如圖1所示。
將連續兩幀IVUS圖像與可變尺度的二維高斯函數G(x,y,σ)卷積運算得到尺度空間L(x,y,σ),對兩個尺度空間劃分模塊后進行匹配。在尺度i計算配準函數的最小化,得到位移場Di,對尺度i插值得到尺度i-1,再對尺度i-1的配準函數最小化,得到新的位移場Di-1,重復相同的操作直到尺度i=0。
其中,i表示第i層尺度空間,i∈(0,...,I)。n表示第n個模塊,n∈(1,,,Ni)。k×l表示每個模塊大小k=2i+2,l=2i。Di表示第i個尺度空間的位移,如式(5):
λi((▽I(x,y)-▽J(x,y))2+Q(▽Di:▽Di)))
(5)

圖2 系統實現流程圖
本系統實現過程需要調整三個關鍵參數:正則化系數,比例閾值和搜索范圍。這些參數通常是根據圖像憑經驗選擇的,本文針對三個參數做了改進,進一步提高方法的性能。
通常在選擇正則項權重λ時采用常數,每個尺度的圖像配準進行相同的正則化導致位移估計過平滑或欠平滑,同時由于原始尺度較大,最小化時易陷入配準函數局部最小值。為此,本文將正則化權重λi進行改進:
對于原始尺度I的正則化權重λi設為零,即λI=0。
對于i
(6)

通常在配準函數最小化過程中采用固定迭代次數,每次迭代都會依次掃描各模塊。然而,不同尺度模塊數量相差較大,導致迭代次數選擇過大時,較大尺度模塊已完成配準,不能立即停止迭代,運行時間較長;迭代次數選擇較小時,劃分精細的尺度未完成配準,位移誤差較大。
為了使迭代均勻化,保證不同尺度圖像既可以完成配準又可以及時停止迭代。定義了一個比例閾值ηi,如式(7),在較大尺度上完成配準后,停止迭代;在模塊劃分精細的尺度上,若未計算的模塊像素數量小于總模塊像素數量的百分之一時,停止迭代。
(7)
其中:選取尺度i=2是因為IVUS圖像尺寸為256×256,尺度空間構建7層,若尺度i選取過大,模塊數量較少使用比例閾值會增加估計誤差;若尺度i選取過小,當模塊數量很大時未能及時停止迭代,閾值效果不明顯。選擇尺度閾值i=2,此尺度以下模塊數量以千量級增加最快。選取百分之一比例保證了加快迭代的同時不影響位移估計準確性。改進的比例閾值法減少了配準過程的迭代次數, 提高了計算速度。
經典多尺度方法在不同尺度下掃描模塊進行配準時,搜索窗口為固定窗口,如式(8)所示:
?i∈{0,...,I},hi=h
(8)
導致估計的位移誤差增加,計算緩慢。
針對IVUS圖像相似性大、位移較小的特點,使用上一級尺度位移來定義下一級尺度的搜索窗口,也就是根據所計算的塊本身來定義搜索窗口,如式(9):
(9)

這種搜索方法使計算位移過程中各模塊區內的點朝向最佳的相鄰點進行計算。如果一個塊的位移估計值在計算過程中沒有被修改,可能意味著這個塊屬于一個同質區域,并且它的估計已經達到一個穩定狀態。改進的搜索方法允許靈活地計算能量函數的最小值,減少計算時間,提高了位移估計準確性。
本文首先利用三類圖像分別是醫院獲得序列IVUS圖像,文獻采用的序列IVUS圖像,人工變形的IVUS圖像,進行了四組實驗,從徑向和角向位移圖驗證結果有效性。
最后從均方根誤差(RMSE)、運行時間、迭代次數等數據方面比較本文方法和改進前方法。
用內外膜分割方法提取IVUS圖像的ROI區域,如圖3所示。利用本文方法和改進前方法進行IVUS圖像配準計算ROI區域位移,結果如圖4和圖5所示。

圖3 ROI區域圖
本文方法的位移結果如圖4所示。

圖4 位移結果圖
改進前方法的位移結果如圖5所示。

圖5 位移結果圖
結果圖顏色軸0以上代表正向位移增加,0以下代表負向位移增加。徑向圖以圖像的中心為圓心,沿徑向方向為正。角向圖以順時針為正,逆時針為負。從結果圖看出改進前和改進后均滿足內膜附近位移較大,外膜附近位移較小的特點,這與實際情況是符合的。但改進前方法的位移結果更粗糙,尤其是角向位移變化很不明顯,誤差較大。本文徑向位移圖和角向位移圖平滑性優于改進前方法,從角向圖明顯血管內外膜位移分布情況。
圖6為文獻[8]采用的連續兩幀IVUS圖像,圖7為本文實驗結果,圖8為改進前的實驗結果:

圖6 連續兩幀IVUS圖像
本文實驗得到的徑向位移和角向位移:

圖7 本文實驗結果
改進前方法的徑向位移和角向位移:

圖8 改進前的實驗結果
利用文獻[8]的IVUS圖像得到的實驗結果圖中同樣看出,改進前的算法內外膜位移變化很小,基本無變化,與實際不太符合。本文位移結果與改進前位移結果相比,清晰看出從內膜區域到外膜區域的位移變化,在十二點到六點區域徑向和角向位移較大,六點和十點區域位移較小,因此本文實驗結果優于改進前的結果。
由于血管變形是無法獲取到的,本文使用液化方法對一幀IVUS圖像內膜12點-1點方向施加壓力。變形前IVUS圖像如圖9(a),變形后IVUS圖像如圖9(b)。其中,圖像的內膜變形大于外膜變形。使用本文方法對變形前和變形后IVUS圖像進行配準,得到圖像ROI區域的徑向和角向位移結果,如圖10(a)和(b)。

圖9 IVUS圖像

圖10 徑向位移和角向位移
從結果圖看出,徑向和角向位移均滿足在12點-1點方向的位移最大,且內膜的位移大于外膜的位移,與人工變形情況一致。
選取斑塊較明顯的IVUS圖像進行配準,對得到的徑向位移圖像進行定量分析。圖11為連續兩幀IVUS圖像,在11點到1點區域和3點區域兩處存在鈣化斑塊。圖12為徑向位移圖。表1為徑向位移量化表。

表1 徑向位移量化表

圖11 平共處五項原則

圖12 向位移
從表格中能夠看出,CD和HG間位移最大,AB和EF間位移最小。與圖11的IVUS圖像比對,AB和EF區域正是IVUS圖像存在鈣化斑塊區域,而CD和GH區域鈣化
斑塊較少。生物力學研究表明[15],鈣化斑塊彈性模量較大,不易變形;纖維和脂質斑塊模量較小,易變性,證明本文得到的位移與生物力學是一致的。
最后,利用RMSE、運行時間、迭代次數等參數比較本文方法和改進前方法。
RMSE參數定義:
(10)
配準實驗進行10次,將10次結果的平均值作為最終結果。對比結果如表2:

表2 本文方法與改進前方法參數比較
從表2看出,本文方法的RMSE,運行時間和迭代次數均低于改進前方法。實驗結果表示,本文方法配準精度提高約25%,運行時間減少20.8 s。
基于IVUS圖像對血管壁位移進行精確估計是一個復雜且有挑戰性的課題,在計算局部組織位移后導出軸向應變就可以獲得組織彈性的量度。通過最小化能量方程來估計位移場,施加灰度平滑約束和光流梯度約束,這樣結合了局部和全局信息的優勢,得到精確并且低噪聲的位移場。尋找全局能量的最小值是非常耗時的,為了加速收斂,使用了從粗到細的多尺度最小化。提出的正則權重,避免了位移估計過平滑;改進的比例閾值,實現了不同尺度迭代次數均勻化,加快了計算速度;結合上一級尺度位移估計值的窗口方法,朝向最佳方向進行搜索配準,提高了位移估計的準確性。
該方法的局限性在于,當血管的角向位移較小或硬化斑塊軸向偏大時,前后兩幀圖像包含較少的角向位移信息,導致計算角向位移存在偏差,這一問題的解決方式是盡量降低導管的回拉速度增加軸向采樣率。
總之,實驗結果表明提出的非剛性圖像配準方法在IVUS圖像位移估計方面的有效性。雖然結果是初步的,但這種方法所得位移圖效果較好,值得進一步深入研究。