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

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

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