基于IVUS圖像配準(zhǔn)的動(dòng)脈壁位移估計(jì)

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

0 引言

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

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

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

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

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

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

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

E=Sim+λReg

(1)

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

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

(2)

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

Reg采用灰度平滑約束和光流梯度約束，用來提高位移估計(jì)的魯棒性，如式(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表示光流梯度,:符號(hào)表示矩陣內(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 多尺度迭代算法及其改進(jìn)

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

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

將連續(xù)兩幀IVUS圖像與可變尺度的二維高斯函數(shù)G(x,y,σ)卷積運(yùn)算得到尺度空間L(x,y,σ)，對(duì)兩個(gè)尺度空間劃分模塊后進(jìn)行匹配。在尺度i計(jì)算配準(zhǔn)函數(shù)的最小化，得到位移場(chǎng)Di，對(duì)尺度i插值得到尺度i-1，再對(duì)尺度i-1的配準(zhǔn)函數(shù)最小化，得到新的位移場(chǎng)Di-1，重復(fù)相同的操作直到尺度i=0。

其中，i表示第i層尺度空間，i∈(0,...,I)。n表示第n個(gè)模塊，n∈(1,,,Ni)。k×l表示每個(gè)模塊大小k=2i+2，l=2i。Di表示第i個(gè)尺度空間的位移，如式(5):

λi((▽I(x,y)-▽J(x,y))2+Q(▽Di:▽Di)))

(5)

圖2 系統(tǒng)實(shí)現(xiàn)流程圖

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

2.1 改進(jìn)的正則項(xiàng)權(quán)重

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

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

對(duì)于i

(6)

2.2 改進(jìn)的比例閾值

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

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

(7)

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

2.3 改進(jìn)的搜索窗口

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

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

(8)

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

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

(9)

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

3 實(shí)驗(yàn)結(jié)果及分析

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

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

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

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

圖3 ROI區(qū)域圖

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

圖4 位移結(jié)果圖

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

圖5 位移結(jié)果圖

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

3.2 基于文獻(xiàn)[8]中IVUS圖像的實(shí)驗(yàn)

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

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

本文實(shí)驗(yàn)得到的徑向位移和角向位移：

圖7 本文實(shí)驗(yàn)結(jié)果

改進(jìn)前方法的徑向位移和角向位移：

圖8 改進(jìn)前的實(shí)驗(yàn)結(jié)果

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

3.3 基于人工變形IVUS圖像的實(shí)驗(yàn)

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

圖9 IVUS圖像

圖10 徑向位移和角向位移

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

3.4 量化位移實(shí)驗(yàn)

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

表1 徑向位移量化表

圖11 平共處五項(xiàng)原則

圖12 向位移

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

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

3.5 參數(shù)比較

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

RMSE參數(shù)定義：

(10)

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

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

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

4結(jié)論

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

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

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