基于多層螺旋計算機(jī)層析數(shù)據(jù)的冠脈三維運(yùn)動估計

基于多層螺旋計算機(jī)層析數(shù)據(jù)的冠脈三維運(yùn)動估計

1血管模型的建立和動力學(xué)冠狀動脈粥樣硬化的估計是心臟疾病研究領(lǐng)域的一個熱點問題之一。x線成像圖像用于跟蹤血管，并分析運(yùn)動的變化。由于X射線造影成像的本質(zhì)是透視投影,不同深度的空間血管段在投影圖像上具有不同的放大率,因此僅從二維圖像估計運(yùn)動是不準(zhǔn)確的。孫正等人提出將三維重建和運(yùn)動估計結(jié)合起來進(jìn)行,由兩幅不同角度的單面二維造影圖像重建出三維血管,之后在三維空間內(nèi)對血管運(yùn)動進(jìn)行分析,但算法的結(jié)果在很大程度上取決于三維重建的精度,而在二維圖像的獲取過程中,平臺的移動不可避免,需要通過優(yōu)化變換矩陣減小重建誤差。近年,多層螺旋計算機(jī)層析(MSCT)成像技術(shù)發(fā)展迅速,特別在冠心病診治中,MSCT冠脈造影(CTA)作為一種無創(chuàng)、三維成像的檢查方式,在臨床上應(yīng)用日益廣泛。不斷提高的時空分辨率,使得利用CTA數(shù)據(jù)分析運(yùn)動的血管成為可能。目前,國內(nèi)外研究多集中在利用CTA數(shù)據(jù)構(gòu)建三維靜態(tài)血管模型進(jìn)行狹窄段分析。本文著眼于動態(tài)模型的建立及運(yùn)動場的估計。對4D-CTA數(shù)據(jù)進(jìn)行分析,利用局部血管增強(qiáng)及自適應(yīng)區(qū)域生長方法提取出各個期相的冠脈樹并進(jìn)行細(xì)化,對骨架點進(jìn)行基于連貫點漂移(CPD)算法的配準(zhǔn),從而估計運(yùn)動場。通過對模擬數(shù)據(jù)和真實全期相CTA數(shù)據(jù)的處理,驗證本文提出方法的精度,并定量描述主要分支的位移幅度與速率。2血管預(yù)分割和溶栓設(shè)置本文在前期工作的基礎(chǔ)上,應(yīng)用基于Hessian矩陣的局部血管增強(qiáng)與自適應(yīng)閾值區(qū)域生長相結(jié)合的方法進(jìn)行血管分割。由于造影劑的作用,CTA數(shù)據(jù)中的血管的灰度通常大于周圍組織,表現(xiàn)為較為明亮的細(xì)長管狀體。對體素數(shù)據(jù)進(jìn)行差分計算,獲得其二階導(dǎo)數(shù)矩陣,即Hessian矩陣。對于血管區(qū)域內(nèi)的點,其Hessian矩陣特征值特點為λ1≈0,λ2<0,λ3<0且|λ2|≈|λ3|?|λ1|,其中λ1對應(yīng)的特征方向為血管方向,通過分析矩陣的特征值,區(qū)分血管和非血管區(qū)域內(nèi)的體素,對數(shù)據(jù)進(jìn)行濾波,從而達(dá)到增強(qiáng)血管、抑制其他組織的目的。為更好地適應(yīng)血管灰度的不均勻性,將數(shù)據(jù)空間劃分為一系列的感興趣區(qū)域(ROI),在每一ROI中先構(gòu)造血管增強(qiáng)濾波函數(shù)進(jìn)行預(yù)處理,再進(jìn)行血管分割。首先指定初始種子點,選擇一個特定大小的ROI,計算該區(qū)域內(nèi)Hessian矩陣的特征值,利用特征值構(gòu)造血管增強(qiáng)濾波函數(shù),抑制非血管結(jié)構(gòu)。對增強(qiáng)后的數(shù)據(jù)進(jìn)行預(yù)分割,從初始值開始,逐次降低閾值,統(tǒng)計分割出的區(qū)域體積,當(dāng)閾值低到一定程度時,分割出的結(jié)果會泄漏到血管周圍區(qū)域,引起一個比較大的體積突變,則上一次的閾值為最佳閾值。將最佳閾值作為最終閾值進(jìn)行區(qū)域生長,利用質(zhì)心法計算分割出的血管在ROI立方體各個截面上的質(zhì)心,作為新的種子點,生成新的ROI,重復(fù)分割過程。每一次分割的局部結(jié)果及時進(jìn)行全局融合,最終得到完整的冠狀動脈樹。ROI的大小由多次實驗確定,經(jīng)驗值為邊長為31voxel的立方體。分割后的冠脈采用基于數(shù)學(xué)形態(tài)學(xué)的三維細(xì)化算法,從三維空間6個方向上進(jìn)行迭代判斷,由外向內(nèi)刪除“簡化點”,保留端點、分支點等特征點,直至得到26連通的單體素曲線,經(jīng)過剪枝等后處理,消除細(xì)小毛刺,獲得保持原血管樹拓?fù)浣Y(jié)構(gòu)及連通特性的骨架線。對各個期相的數(shù)據(jù)均進(jìn)行同樣的冠脈分割及細(xì)化處理步驟,得到同一冠脈分支不同時刻的骨架點序列,記為Xt={xt,i}(i=1,2,……nt),t為表征不同時刻的序號,nt為該時刻分支上點的個數(shù)。3高斯模型權(quán)重的確定對于運(yùn)動場的估計,一般轉(zhuǎn)化為點集之間的配準(zhǔn)問題進(jìn)行求解。假設(shè)點集Xt與Xt+1分別為相鄰兩個時刻某分支上的骨架點,為方便,Xt記為Y,Xt+1記為X。對兩點集配準(zhǔn),就是確定兩個點集之間的對應(yīng)點對{xi,yj}及空間位置變換關(guān)系T,使得xi=T(yj)。圖像配準(zhǔn)算法很多,具體到冠脈運(yùn)動,因為軌跡復(fù)雜,且在運(yùn)動過程中發(fā)生形變,所以T為非剛體變換。較為常見的非剛體配準(zhǔn)方法有自由形變模型(FFD)、魯棒點配準(zhǔn)(RPM)和CPD算法等,其中CPD算法將運(yùn)動連貫理論(MCT)應(yīng)用于配準(zhǔn)過程,保證了相鄰點之間變換的平滑性,在非剛體點配準(zhǔn)中取得了較為滿意的效果。假設(shè)X與Y分別為數(shù)據(jù)點集與模型點集,記為XN×D=(x1,x2,…xN),YM×D=(y1,y2,…yM),維度D=3。記Y變換后的點集為T(Y,θ),θ為變換參數(shù)。若將X與T(Y,θ)視為高斯混合模型的點分布與相應(yīng)的模型中心,則X的分布服從p(xn)=Μ+1∑i=1Ρ(m)p(xn|m)=w1Ν+(1-w)Μ∑i=1Ρ(m)1(2πσ2)D/2exp(-∥xn-Τ(ym?θ)∥22σ2)(1)p(xn)=∑i=1M+1P(m)p(xn|m)=w1N+(1?w)∑i=1MP(m)1(2πσ2)D/2exp(?∥xn?T(ym?θ)∥22σ2)(1)其中:P(m)為各個高斯模型的權(quán)重,簡化起見,規(guī)定P(m)=1/M;w為無匹配點在點集中所占的比例;σ2為高斯模型的協(xié)方差。xn與ym間的對應(yīng)關(guān)系由后驗概率計算,即P(m|xn)=P(m)p(xn|m)/p(xn)(2)式(1)中,未知參數(shù)為θ、σ2。構(gòu)造目標(biāo)函數(shù),當(dāng)函數(shù)具有最小值時的參數(shù)取值為最優(yōu)解,即E(θ?σ2)=-Ν∑n=1logΜ+1∑m=1Ρ(m)p(xn|m)(3)E(θ?σ2)=?∑n=1Nlog∑m=1M+1P(m)p(xn|m)(3)改寫式(3),利用期望最大化方法(EM)算法求解Q(θ?σ2)=12σ2Μ?Ν∑m?n=1Ρ(m|xn)∥xn-(ym+Τ(ym))∥2+ΝΡD2logσ2(4)Q(θ?σ2)=12σ2∑m?n=1M?NP(m|xn)∥xn?(ym+T(ym))∥2+NPD2logσ2(4)其中,Νp=Ν∑n=1Μ∑m=1Ρ(m|xn)Np=∑n=1N∑m=1MP(m|xn)對于剛體變換,式(4)中T可進(jìn)一步分解為平移、旋轉(zhuǎn)和縮放等矩陣;而對于非剛體變換則需要選擇更為復(fù)雜的函數(shù)形式v,對形變進(jìn)行滿足平滑性和準(zhǔn)確性的逼近。對式(4)加入正則化項,進(jìn)行平滑性約束,有Q(v?σ2)=12σ2Μ?Ν∑m?n=1Ρold(m|xn)∥xn-(ym+v(ym))∥2+ΝΡD2logσ2+λ2∥Ρv∥(5)Q(v?σ2)=12σ2∑m?n=1M?NPold(m|xn)∥xn?(ym+v(ym))∥2+NPD2logσ2+λ2∥Pv∥(5)根據(jù)MCT理論,相鄰點的運(yùn)動是連貫的,即運(yùn)動場v是平滑的,若將v用參數(shù)為β的高斯核的線性組合表示,則正好滿足這種平滑性要求,即v(yi)=Μ∑m=1wmexp(-12∥yi-ymβ∥)=1σ2λΝ∑n=1Ρ(mxn)(xn-(ym+v(ym))exp(-12)∥yi-ymβ∥)(6)v(yi)=∑m=1Mwmexp(?12∥yi?ymβ∥)=1σ2λ∑n=1NP(mxn)(xn?(ym+v(ym))exp(?12)∥yi?ymβ∥)(6)求解最小二乘問題及變分方程,得到v及σ2的最優(yōu)解;再根據(jù)式(2),計算點集間的對應(yīng)關(guān)系矩陣,找到使得后驗概率最大的對應(yīng)點對,從而得到點集Y到X的形變場。4結(jié)果與分析4.1oxel1333實驗數(shù)據(jù)為64層螺旋CT機(jī)掃描圖像,分為11個期相,每個期相的數(shù)據(jù)均為512voxel×512voxel×215voxel,層間分辨率為0.625mm,層內(nèi)分辨率為0.3899mm×0.3899mm。圖1是單期相分割出的左冠脈及細(xì)化后的骨架。將多期相的左、右冠脈都分割出來,并進(jìn)行細(xì)化,骨架線如圖2所示。4.2位移估計誤差為了定量評價配準(zhǔn)算法精度,對某一時刻的回旋支LCX骨架構(gòu)造已知參數(shù)的運(yùn)動場,得到形變后的曲線,然后對這兩條曲線進(jìn)行配準(zhǔn)進(jìn)而估計運(yùn)動場,比較估計出的運(yùn)動場與已知運(yùn)動場的參數(shù),評價算法的誤差。心臟的運(yùn)動主要為非剛體的膨脹和收縮,可近似用“旋轉(zhuǎn)+平移+放大/收縮+扭轉(zhuǎn)”模擬。假設(shè)Rz、Ry和Rx分別表示繞三維空間z、y、x3個軸的旋轉(zhuǎn)矩陣,E表示縮放矩陣,S表示扭轉(zhuǎn)變形矩陣,T為平移矩陣,則某點y(y1,y2,y3)運(yùn)動后的坐標(biāo)x(x1,x2,x3)可表示為[x1x2x3]=RzRyRxES[y1y2y3]+Τ=[cosγsinγ0-sinγcosγ0001][cosβ0-sinβ010sinβ0cosβ][1000cosαsinα0-sinαcosα][Ex000Ey000Ez][1Κxx-Κyy-Κxx1ΚzzΚyy-Κzz][y1y2y3]+[ΤxΤyΤz](7)式(7)中,各項參數(shù)的取值見表1。膨脹、收縮時構(gòu)造出的運(yùn)動場及估計出的運(yùn)動場如圖3、圖4所示。直觀上,圖3、圖4估計出的運(yùn)動場都是比較平滑和均勻的,且和構(gòu)造的運(yùn)動場比較相似。對兩種情況下位移情況進(jìn)行計算,其中Drmax、Drmin、Dravex、Demax、Demin和Deavex分別表示構(gòu)造的運(yùn)動場與估計出的運(yùn)動場的最大位移、最小位移和平均位移,見表2。以最大位移為例,定義最大位移估計誤差為Emax=Demax-DrmaxDrmax×100％(8)其他兩種位移誤差定義與式(8)類似,配準(zhǔn)誤差以均方根誤差(RMSE)衡量。兩種情況下的配準(zhǔn)RMSE分別為0.4744mm和0.3512mm,小于1voxel,位移估計誤差小于1%。由計算結(jié)果可知,估計出的運(yùn)動場與構(gòu)造的運(yùn)動場差別很小,采用CPD進(jìn)行運(yùn)動估計,具有較高的精度。4.3降壓降位計算時各各自為基的平均移動距離、平均位移幅度及速率對全期相的左右冠脈主要分支前降支LAD、回旋支LCX、對角支D及右冠脈RCA分別進(jìn)行配準(zhǔn),并計算其在各個時刻的3D位移幅度和速率。位移幅度定義為其他時刻相對于R-R間期0%時各分支上點的平均移動距離;速度定義為Vi=Di/t(i=1,2,……10),其中Di是該期相到下一期相的平均位移幅度,t為期相間的時間間隔,以心率60bpm計算,則t=0.1s。計算結(jié)果見表3、表4。對表3、4中的數(shù)據(jù)進(jìn)行三次樣條插值,得到整個R-R間期內(nèi)任意時刻的各分支的3D平均位移幅度及速率,如圖5、圖6所示。從圖可以看出,右冠脈相對于左冠脈運(yùn)動更為明顯,具有更大的位移幅度及更高的速率,同時左、右冠脈均有3個速率較高的時刻,對應(yīng)心臟運(yùn)動中的收縮、舒張期充盈及心房收