(北京工業(yè)大學(xué) 信息學(xué)部，北京 100124)

0 引言

隨著年齡的增長，心血管事件如冠狀動脈疾病或動脈粥樣硬化，常常引起生物組織硬度變化和動脈粥樣硬化斑塊的產(chǎn)生，最終導(dǎo)致中風(fēng)或心肌梗死。最初，將安裝超聲儀的導(dǎo)管進入血管內(nèi)獲得IVUS成像，來識別動脈壁形態(tài)，纖維斑塊位置和厚度以及內(nèi)腔面積，改善心血管疾病(如動脈粥樣硬化)的診斷。因此，其臨床應(yīng)用局限于輔助性的超聲波檢查工具或精湛的血管內(nèi)手術(shù)。

然而，彈性成像是一種非常有前景的方法來表征血管壁的力學(xué)性能。分析血管彈性對預(yù)防斑塊硬化、易損斑塊的早期發(fā)現(xiàn)以及指導(dǎo)動脈治療手術(shù)具有十分重要的意義[1]。生物力學(xué)研究表明[3]，動脈粥樣硬化斑塊某個位置的應(yīng)力超過局部閾值時易發(fā)生破裂。但是，如何準(zhǔn)確評價血管斑塊彈性及預(yù)測斑塊破裂風(fēng)險的方法仍然缺乏[4]。

非剛性圖像配準(zhǔn)一直是醫(yī)學(xué)圖像處理領(lǐng)域研究的熱點[5]，可以用來計算血管壁變形，研究斑塊破裂風(fēng)險，從而為預(yù)防和治療心血管疾病提供依據(jù)[6-7]。LiangYun等人[7]使用有限元網(wǎng)格和三次Bspline基函數(shù)作為配準(zhǔn)模型，采用列文伯格-馬夸爾特(LM)法方法最小化配準(zhǔn)函數(shù)。然而，三次Bsplines的離散初始化過分平滑應(yīng)變值，使測量的應(yīng)變不準(zhǔn)確。Richards等人[9-10]使用半范數(shù)正則化配準(zhǔn)模型，采用牛頓下降法最小化配準(zhǔn)函數(shù)。但位移較大時，變形效果不太真實。Tang Z[11]等人在配準(zhǔn)能模型引入平衡系數(shù)，一定程度緩解了配準(zhǔn)精度較低問題，但計算速度較慢。

由于處理IVUS圖像本身不足以解決由于高噪聲造成的位移模糊問題，因此使用正則化方法為整合這些約束信息提供了必要的理論框架。由于血管內(nèi)外膜運動的復(fù)雜性和空間異質(zhì)性，全局參數(shù)正則化模型是不合適的。而光流方法是一個更好的模型，它允許估計局部變形。因此，為了進一度提高位移估計的配準(zhǔn)精度和計算速度，本文提出一種結(jié)合光流的非剛性配準(zhǔn)模型，使用改進的多尺度方法最小化配準(zhǔn)函數(shù)，得到了較好的位移結(jié)果。

1 改進的圖像配準(zhǔn)模型

非剛性圖像配準(zhǔn)被認為是一種優(yōu)化問題，其目標(biāo)是最小化配準(zhǔn)函數(shù)。 在本文中，所提出的配準(zhǔn)函數(shù)包括表示前后兩幀IVUS圖像的相似性項，以及結(jié)合組織位移平滑約束增加位移魯棒性的加權(quán)項。

經(jīng)典的配準(zhǔn)模型[12]，如式(1)：

E=Sim+λReg

(1)

E表示圖像配準(zhǔn)函數(shù)，Sim表示圖像配準(zhǔn)項，用來評價前后兩幀圖像的相關(guān)性。Reg表示正則項，用來平滑位移函數(shù)。λ表示正則項比重。

其中，Sim采用平方和差(SSD)準(zhǔn)則[12]，，如式(2)：

(2)

I(x,y)代表前一幀IVUS圖像,J(x,y)代表后一幀IVUS圖像。D表示圖像I和圖像J的變換矩陣，即待求解的位移場。Ω表示IVUS圖像中計算位移的區(qū)域。

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表示光流梯度,:符號表示矩陣內(nèi)積，(▽D:▽D)表示光流梯度約束。Q表示空間權(quán)重函數(shù)，用于調(diào)節(jié)光流梯度約束，本文以血管內(nèi)膜距血管中心的徑向距離R(x,y)表示即Q=R4。

將式(2)和式(3)代入式(1)中得到最終配準(zhǔn)模型，如式(4)：

(4)

2 多尺度迭代算法及其改進

最小化能量函數(shù)是一個經(jīng)典的優(yōu)化問題,能量函數(shù)是非凸的，可能包含幾個局部最小值。其次，最小化計算時間是一個重要的考慮因素。如模擬退火等隨機算法在理論上向能量函數(shù)的全局最小值收斂。然而，隨機松弛算法需要大量的迭代非常緩慢；迭代條件模式(ICM)算法等確定性算法會更快地收斂，但可以能陷入能量函數(shù)的局部極小值中。為了避免選擇局部最小值，同時減少迭代時間，多尺度方法通過在更精細的配置空間內(nèi)搜索提供了這種可能性。

多尺度迭代方法[13-14]如圖1所示。

將連續(xù)兩幀IVUS圖像與可變尺度的二維高斯函數(shù)G(x,y,σ)卷積運算得到尺度空間L(x,y,σ)，對兩個尺度空間劃分模塊后進行匹配。在尺度i計算配準(zhǔn)函數(shù)的最小化，得到位移場Di，對尺度i插值得到尺度i-1，再對尺度i-1的配準(zhǔn)函數(shù)最小化，得到新的位移場Di-1，重復(fù)相同的操作直到尺度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 系統(tǒng)實現(xiàn)流程圖

本系統(tǒng)實現(xiàn)過程需要調(diào)整三個關(guān)鍵參數(shù)：正則化系數(shù)，比例閾值和搜索范圍。這些參數(shù)通常是根據(jù)圖像憑經(jīng)驗選擇的，本文針對三個參數(shù)做了改進，進一步提高方法的性能。

2.1 改進的正則項權(quán)重

通常在選擇正則項權(quán)重λ時采用常數(shù)，每個尺度的圖像配準(zhǔn)進行相同的正則化導(dǎo)致位移估計過平滑或欠平滑，同時由于原始尺度較大，最小化時易陷入配準(zhǔn)函數(shù)局部最小值。為此，本文將正則化權(quán)重λi進行改進：

對于原始尺度I的正則化權(quán)重λi設(shè)為零，即λI=0。

對于i

(6)

2.2 改進的比例閾值

通常在配準(zhǔn)函數(shù)最小化過程中采用固定迭代次數(shù)，每次迭代都會依次掃描各模塊。然而，不同尺度模塊數(shù)量相差較大，導(dǎo)致迭代次數(shù)選擇過大時，較大尺度模塊已完成配準(zhǔn)，不能立即停止迭代，運行時間較長；迭代次數(shù)選擇較小時，劃分精細的尺度未完成配準(zhǔn)，位移誤差較大。

為了使迭代均勻化，保證不同尺度圖像既可以完成配準(zhǔn)又可以及時停止迭代。定義了一個比例閾值ηi，如式(7)，在較大尺度上完成配準(zhǔn)后，停止迭代；在模塊劃分精細的尺度上，若未計算的模塊像素數(shù)量小于總模塊像素數(shù)量的百分之一時，停止迭代。

(7)

其中：選取尺度i=2是因為IVUS圖像尺寸為256×256，尺度空間構(gòu)建7層，若尺度i選取過大，模塊數(shù)量較少使用比例閾值會增加估計誤差；若尺度i選取過小，當(dāng)模塊數(shù)量很大時未能及時停止迭代，閾值效果不明顯。選擇尺度閾值i=2，此尺度以下模塊數(shù)量以千量級增加最快。選取百分之一比例保證了加快迭代的同時不影響位移估計準(zhǔn)確性。改進的比例閾值法減少了配準(zhǔn)過程的迭代次數(shù)， 提高了計算速度。

2.3 改進的搜索窗口

經(jīng)典多尺度方法在不同尺度下掃描模塊進行配準(zhǔn)時，搜索窗口為固定窗口，如式(8)所示:

?i∈{0,...,I},hi=h

(8)

導(dǎo)致估計的位移誤差增加，計算緩慢。

針對IVUS圖像相似性大、位移較小的特點，使用上一級尺度位移來定義下一級尺度的搜索窗口，也就是根據(jù)所計算的塊本身來定義搜索窗口，如式(9):

(9)

這種搜索方法使計算位移過程中各模塊區(qū)內(nèi)的點朝向最佳的相鄰點進行計算。如果一個塊的位移估計值在計算過程中沒有被修改，可能意味著這個塊屬于一個同質(zhì)區(qū)域，并且它的估計已經(jīng)達到一個穩(wěn)定狀態(tài)。改進的搜索方法允許靈活地計算能量函數(shù)的最小值，減少計算時間，提高了位移估計準(zhǔn)確性。

3 實驗結(jié)果及分析

本文首先利用三類圖像分別是醫(yī)院獲得序列IVUS圖像，文獻采用的序列IVUS圖像，人工變形的IVUS圖像，進行了四組實驗，從徑向和角向位移圖驗證結(jié)果有效性。

最后從均方根誤差(RMSE)、運行時間、迭代次數(shù)等數(shù)據(jù)方面比較本文方法和改進前方法。

3.1 基于醫(yī)院的連續(xù)兩幀IVUS圖像配準(zhǔn)實驗

用內(nèi)外膜分割方法提取IVUS圖像的ROI區(qū)域，如圖3所示。利用本文方法和改進前方法進行IVUS圖像配準(zhǔn)計算ROI區(qū)域位移，結(jié)果如圖4和圖5所示。

圖3 ROI區(qū)域圖

本文方法的位移結(jié)果如圖4所示。

圖4 位移結(jié)果圖

改進前方法的位移結(jié)果如圖5所示。

圖5 位移結(jié)果圖

結(jié)果圖顏色軸0以上代表正向位移增加，0以下代表負向位移增加。徑向圖以圖像的中心為圓心，沿徑向方向為正。角向圖以順時針為正，逆時針為負。從結(jié)果圖看出改進前和改進后均滿足內(nèi)膜附近位移較大，外膜附近位移較小的特點，這與實際情況是符合的。但改進前方法的位移結(jié)果更粗糙，尤其是角向位移變化很不明顯，誤差較大。本文徑向位移圖和角向位移圖平滑性優(yōu)于改進前方法，從角向圖明顯血管內(nèi)外膜位移分布情況。

3.2 基于文獻[8]中IVUS圖像的實驗

圖6為文獻[8]采用的連續(xù)兩幀IVUS圖像，圖7為本文實驗結(jié)果，圖8為改進前的實驗結(jié)果：

圖6 連續(xù)兩幀IVUS圖像

本文實驗得到的徑向位移和角向位移：

圖7 本文實驗結(jié)果

改進前方法的徑向位移和角向位移：

圖8 改進前的實驗結(jié)果

利用文獻[8]的IVUS圖像得到的實驗結(jié)果圖中同樣看出，改進前的算法內(nèi)外膜位移變化很小，基本無變化，與實際不太符合。本文位移結(jié)果與改進前位移結(jié)果相比，清晰看出從內(nèi)膜區(qū)域到外膜區(qū)域的位移變化，在十二點到六點區(qū)域徑向和角向位移較大，六點和十點區(qū)域位移較小，因此本文實驗結(jié)果優(yōu)于改進前的結(jié)果。

3.3 基于人工變形IVUS圖像的實驗

由于血管變形是無法獲取到的，本文使用液化方法對一幀IVUS圖像內(nèi)膜12點-1點方向施加壓力。變形前IVUS圖像如圖9(a)，變形后IVUS圖像如圖9(b)。其中，圖像的內(nèi)膜變形大于外膜變形。使用本文方法對變形前和變形后IVUS圖像進行配準(zhǔn)，得到圖像ROI區(qū)域的徑向和角向位移結(jié)果，如圖10(a)和(b)。

圖9 IVUS圖像

圖10 徑向位移和角向位移

從結(jié)果圖看出，徑向和角向位移均滿足在12點-1點方向的位移最大，且內(nèi)膜的位移大于外膜的位移，與人工變形情況一致。

3.4 量化位移實驗

選取斑塊較明顯的IVUS圖像進行配準(zhǔn)，對得到的徑向位移圖像進行定量分析。圖11為連續(xù)兩幀IVUS圖像，在11點到1點區(qū)域和3點區(qū)域兩處存在鈣化斑塊。圖12為徑向位移圖。表1為徑向位移量化表。

表1 徑向位移量化表

圖11 平共處五項原則

圖12 向位移

從表格中能夠看出，CD和HG間位移最大，AB和EF間位移最小。與圖11的IVUS圖像比對，AB和EF區(qū)域正是IVUS圖像存在鈣化斑塊區(qū)域，而CD和GH區(qū)域鈣化

斑塊較少。生物力學(xué)研究表明[15]，鈣化斑塊彈性模量較大，不易變形；纖維和脂質(zhì)斑塊模量較小，易變性，證明本文得到的位移與生物力學(xué)是一致的。

3.5 參數(shù)比較

最后，利用RMSE、運行時間、迭代次數(shù)等參數(shù)比較本文方法和改進前方法。

RMSE參數(shù)定義：

(10)

配準(zhǔn)實驗進行10次，將10次結(jié)果的平均值作為最終結(jié)果。對比結(jié)果如表2：

表2 本文方法與改進前方法參數(shù)比較

從表2看出，本文方法的RMSE，運行時間和迭代次數(shù)均低于改進前方法。實驗結(jié)果表示，本文方法配準(zhǔn)精度提高約25%，運行時間減少20.8 s。

4結(jié)論

基于IVUS圖像對血管壁位移進行精確估計是一個復(fù)雜且有挑戰(zhàn)性的課題，在計算局部組織位移后導(dǎo)出軸向應(yīng)變就可以獲得組織彈性的量度。通過最小化能量方程來估計位移場，施加灰度平滑約束和光流梯度約束，這樣結(jié)合了局部和全局信息的優(yōu)勢，得到精確并且低噪聲的位移場。尋找全局能量的最小值是非常耗時的，為了加速收斂，使用了從粗到細的多尺度最小化。提出的正則權(quán)重，避免了位移估計過平滑；改進的比例閾值，實現(xiàn)了不同尺度迭代次數(shù)均勻化，加快了計算速度；結(jié)合上一級尺度位移估計值的窗口方法，朝向最佳方向進行搜索配準(zhǔn)，提高了位移估計的準(zhǔn)確性。

該方法的局限性在于，當(dāng)血管的角向位移較小或硬化斑塊軸向偏大時，前后兩幀圖像包含較少的角向位移信息，導(dǎo)致計算角向位移存在偏差，這一問題的解決方式是盡量降低導(dǎo)管的回拉速度增加軸向采樣率。

總之，實驗結(jié)果表明提出的非剛性圖像配準(zhǔn)方法在IVUS圖像位移估計方面的有效性。雖然結(jié)果是初步的，但這種方法所得位移圖效果較好，值得進一步深入研究。