專利名稱：一種冠狀動脈血管軸線的四維重建方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種根據(jù)覆蓋一個或多個心動周期的數(shù)字X射線冠狀動脈造影圖像序列對血管軸線進行四維(三維+時間)重建的方法，屬醫(yī)學(xué)檢測技術(shù)領(lǐng)域。

背景技術(shù)：
X射線冠狀動脈造影(X-ray coronary angiography，CAG)是目前臨床廣泛采用的診斷和治療冠心病的影像手段。CAG診斷的最大特點在于可視性，即通過靜止或動態(tài)觀察造影劑的充盈和消失情況來判斷心血管解剖學(xué)形態(tài)異常的部位、性質(zhì)和程度。X射線造影成像是將空間血管結(jié)構(gòu)重疊到二維成像平面上，即投影成像，從而丟失了臨床診斷中所需的大部分三維空間信息。傳統(tǒng)的診斷過程中，醫(yī)生需要利用解剖、病理等專門知識和臨床經(jīng)驗，根據(jù)血管在多個方向的投影想象其三維形態(tài)，從而做出評價。分析結(jié)果的準確程度密切依賴于醫(yī)生的臨床經(jīng)驗和專業(yè)知識，客觀性和可重復(fù)性差。同時，冠狀動脈附著在心外膜表面上，隨著心臟有節(jié)律地收縮和舒張，因而在心動周期的不同時刻，血管的形態(tài)結(jié)構(gòu)和空間位置會發(fā)生很大變化。單純依靠一個時刻的二維投影或三維重建結(jié)果很難得到對血管形態(tài)結(jié)構(gòu)的準確描述。
在一個或多個心動周期中拍攝的X射線冠狀動脈造影圖像序列記錄了心臟搏動時的血管變形和心臟供血功能等動態(tài)信息。采用數(shù)字圖像處理的方法，可以從不同時刻、不同角度拍攝的圖像序列中客觀、定量地提取隱含的空間和時間信息，并以恰當(dāng)?shù)姆绞奖磉_出來。這些信息對于冠心病的臨床診治具有重要的輔助作用。同時冠脈的四維(三維+時間)重建也是造影圖像與其它冠脈成像技術(shù)(如血管內(nèi)超聲(IVUS)、血管內(nèi)OCT等)的融合、心動周期中的冠脈運動估計等的關(guān)鍵步驟。
為了實現(xiàn)對冠脈在整個心動周期中的四維重建，現(xiàn)有的方法多是按照“自底向上”，即2D→3D的順序，逐個時刻進行三維重建。首先從兩個角度的造影圖像中提取出主要血管分支的軸線；之后，進行左右兩個角度間軸線對應(yīng)點的匹配；最后，根據(jù)兩個投影點的坐標(biāo)，計算出各軸線點的三維坐標(biāo)。由于需要首先從原始圖像中提取出單像素、八連通的血管軸線，因此結(jié)果的精度在很大程度上取決于二維提取的準確度。但臨床采集的造影圖像通常有較嚴重的噪聲污染，而且圖像中經(jīng)常會出現(xiàn)血管重疊的現(xiàn)象，這種情況下僅從一個角度的圖像很難得到準確的血管輪廓。目前全自動的提取算法尚不能對何種質(zhì)量的圖像均保成功，仍需操作者的手動參與。同時，兩個角度之間血管軸線投影點的匹配直接影響到重建的精度。對于離散的二維軸線，一般是利用立體成像中的外極線約束逐點匹配，該方法很難達到較高的精度，重建結(jié)果的連續(xù)性不好，計算量也很大。此外，由于需要對序列中每個時刻的圖像均重復(fù)上述步驟，因此要完成整個序列的血管跟蹤，所需的計算成本過高，且需要操作者大量的手動參與，實用價值較低。


發(fā)明內(nèi)容
本發(fā)明的目的是克服現(xiàn)有技術(shù)的不足、提供一種操作簡便、計算量小、工作效率高的冠狀動脈血管軸線的四維重建方法。
本發(fā)明所稱問題是以下述技術(shù)方案實現(xiàn)的 一種冠狀動脈血管軸線的四維重建方法，它是在采集兩個角度的覆蓋一個或多個心動周期的X射線冠狀動脈造影圖像序列后，首先建立X射線血管造影系統(tǒng)在兩個角度的投影模型，并推導(dǎo)出兩幅不同角度的造影圖像之間的幾何變換關(guān)系，再對手動選取的第一時刻圖像中感興趣血管段的采樣點進行三維重建，并以連接各重建點所得折線作為snake模型的初始位置，通過snake變形獲得該段血管第一時刻的三維軸線，而對于圖像序列中的每一后續(xù)時刻，則均以其前一時刻的3D血管軸線作為snake模型的初始位置，通過snake變形獲得三維血管軸線，從而完成整個序列血管軸線的四維重建，具體步驟如下 a、以心電信號為同步信號，采集兩個角度的、覆蓋一個或多個心動周期的單面X射線冠狀動脈造影同步圖像序列，并對原始圖像進行初步的濾波、去噪、畸變校正等處理，增強視覺效果； b、建立X射線血管造影系統(tǒng)在兩個角度的投影模型，推導(dǎo)兩幅不同角度的造影圖像之間的幾何變換關(guān)系 設(shè)點s1和s2為兩次造影過程中X射線源焦點的位置，分別以s1和s2為原點，建立空間坐標(biāo)系s1x1y1z1和s2x2y2z2；坐標(biāo)系U1V1O1和U2V2O2為成像平面坐標(biāo)系；D1和D2分別是s1和s2到各自成像平面的垂直距離，空間血管上的點P在成像平面A和B上的投影點分別為p1(u1，v1)和p2(u2，v2)， 則s1x1y1z1和s2x2y2z2之間的幾何變換關(guān)系為 其中，R3×3為旋轉(zhuǎn)矩陣R＝RY(α2)·RX(β2)·RX(-β1)·RY(α1)

為平移向量 L1和L2為X射線源焦點s1和s2到旋轉(zhuǎn)中心的距離；(α1，β1)和(α2，β2)分別是A和B的造影角度；(x1，y1，z1)和(x2，y2，z2)分別表示空間點P在坐標(biāo)系s1x1y1z1和s2x2y2z2中的坐標(biāo)； 根據(jù)透視投影成像的幾何關(guān)系可知，空間點的三維坐標(biāo)與其二維投影坐標(biāo)之間的關(guān)系為 其中， c、在序列中第一時刻的圖像中，通過人工取點獲得感興趣血管段在左右投影中的近似中心線，用折線表示，并利用外極約束得到兩個角度間對應(yīng)點的匹配； d、對手動選取的采樣點進行三維重建； e、連接各3D點，所得折線作為3D snake的初始位置，通過使預(yù)先設(shè)定的能量函數(shù)最小，snake在空間中變形，最終得到具有最小能量的最優(yōu)位置，就是第一時刻的3D血管軸線； f、后續(xù)時刻血管軸線的四維重建 對于序列中的后續(xù)時刻，以前一時刻的3D血管軸線作為當(dāng)前時刻snake的初始位置，通過snake變形，完成整個序列中血管軸線的四維重建。
上述冠狀動脈血管軸線的四維重建方法，所述X射線冠狀動脈造影圖像序列的兩個采集角度之間夾角的取值范圍為60°～120°。
上述冠狀動脈血管軸線的四維重建方法，在第一時刻圖像的感興趣血管段中選取的三維重建采樣點包括起點、終點和3～6個中間點。
本發(fā)明按照“自頂向下”的方式，采用3D snake模型技術(shù)，通過使能量函數(shù)最小，表示血管軸線的曲線直接在三維空間中變形，完成血管軸線的四維重建。該方法將序列圖像中血管的二維中心線提取、三維重建和相鄰時刻間的三維運動跟蹤集成到一個框架中完成，提高了運算精度和速度。與“自底向上”的傳統(tǒng)重建算法相比較，由于曲線直接在三維空間中變形，因此避免了逐點匹配，提高了重建精度。同時操作者的參與被減小到只需在序列的首幀中選擇感興趣血管段上的若干點。本發(fā)明不僅大大減少了由操作者參與所引入的不確定性和誤差，提高了結(jié)果的可重復(fù)性，而且操作簡便，工作效率高。



圖1是本發(fā)明的實施步驟流程圖； 圖2是本發(fā)明的單面造影系統(tǒng)在兩個角度的成像示意圖； 圖3是本發(fā)明的血管造影系統(tǒng)的造影角度示意圖； 圖4是采用本發(fā)明的左冠前降支(LAD)首時刻軸線重建的實施例圖像；其中，圖4(a)是首時刻的造影圖像對和初始的3D snake；圖4(b)是首時刻的造影圖像對、重建出的血管軸線及其在兩個成像平面上的二維投影； 圖5是采用本發(fā)明的左冠前降支(LAD)后續(xù)時刻軸線序列重建的實施例圖像，其中，圖5(a)～(d)分別是序列中的第4、6、8、10個時刻的血管軸線重建結(jié)果及其在兩個成像平面上的二維投影。
圖中各符號為A、B、成像平面；s1、s2、兩次造影過程中X射線源焦點的位置；s1x1y1z1、以s1為原點的空間坐標(biāo)系；s2x2y2z2、以s2為原點的空間坐標(biāo)系；OXYZ、以旋轉(zhuǎn)中心為原點的空間坐標(biāo)系；U1V1O1、成像平面A上的直角坐標(biāo)系；U2V2O2、成像平面B上的直角坐標(biāo)系；D1、s1到成像平面A的垂直距離；D2、s2到成像平面B的垂直距離；P、空間血管上的點；p1、P點在成像平面A上的投影；p2、P點在成像平面B上的投影；u1、p1在坐標(biāo)系U1V1O1內(nèi)的橫坐標(biāo)；v1、p1在坐標(biāo)系U1V1O1內(nèi)的縱坐標(biāo)；u2、p2在坐標(biāo)系U2V2O2內(nèi)的橫坐標(biāo)；v1、p2在坐標(biāo)系U2V2O2內(nèi)的縱坐標(biāo)；K1、K2、外極線；(α1，β1)、成像平面A的造影角度；(α2，β2)、成像平面B的造影角度；L1、X射線源s1到旋轉(zhuǎn)中心的距離；L2、X射線源s2到旋轉(zhuǎn)中心的距離。

具體實施例方式 下面結(jié)合附圖和實例對本發(fā)明作進一步詳細說明。
如附圖1所示，本發(fā)明方法的步驟包括 (1)圖像采集和預(yù)處理 采用C型臂單面X射線血管造影系統(tǒng)獲得兩個角度的、至少覆蓋一個心動周期的冠狀動脈造影圖像序列，要求兩個角度之間的夾角在60°至120°之間。記錄成像系統(tǒng)參數(shù)(造影角度，X射線源到成像平面的距離)。采用同步記錄的心電信號選取心動周期中相同相位的兩個角度的單面圖像對。
對原始圖像進行必要的預(yù)處理，主要包括畸變校正、均衡對比度、去除噪聲和圖像增強等，提高圖像的視覺效果，為后續(xù)處理奠定基礎(chǔ)。
(2)建立單面造影系統(tǒng)在兩個近似正交角度的成像模型 如附圖2所示，點s1和s2表示兩次造影過程中X射線源焦點的位置。分別以s1和s2為原點，建立空間坐標(biāo)系s1x1y1z1和s2x2y2z2；坐標(biāo)系U1V1O1和U2V2O2為成像平面坐標(biāo)系；D1和D2分別是s1和s2到各自成像平面的垂直距離，隨成像面的移動而改變；空間血管上的點P在成像平面A和B上的投影點分別為p1(u1，v1)和p2(u2，v2)。
從坐標(biāo)系s1x1y1z1到s2x2y2z2的變換運動是三維空間中的剛體運動，根據(jù)剛體運動理論，推導(dǎo)出坐標(biāo)系s1x1y1z1和s2x2y2z2之間的幾何變換關(guān)系 其中R3×3為旋轉(zhuǎn)矩陣 R＝RY(α2)·RX(β2)·RX(-β1)·RY(α1) (2)

為平移向量 式中，L1和L2為X射線源s1和s2到旋轉(zhuǎn)中心的距離；(α1，β1)和(α2，β2)分別是圖像A和B的造影角度(如附圖3所示)。距離和角度值均可從造影機上得到。將稱為幾何變換矩陣。
根據(jù)透視投影成像的幾何關(guān)系可知 其中 (x1，y1，z1)和(x2，y2，z2)分別表示空間點P在坐標(biāo)系s1x1y1z1和s2x2y2z2中的坐標(biāo)。因此根據(jù)式(1)、(4)和(5)，由點p1和p2的坐標(biāo)可反算出點P的三維坐標(biāo)。
(3)第一時刻血管軸線的三維重建 在跟蹤過程開始前，需要確定snake模型的初始位置。本發(fā)明采用手動取點的方法，首先從一個角度的圖像中手動選取感興趣血管段上的若干采樣點(一般選取血管段的起點、終點和3～6個中間點即可，具體數(shù)目根據(jù)血管段的長度決定)。然后采用外極約束得到各點在另一個角度圖像上的對應(yīng)點。如附圖2所示，空間點P與X射線源焦點s1和s2構(gòu)成外極平面，該平面與圖像平面A和B相交形成左右兩條外極線K1和K2。根據(jù)外極約束原理，點P在圖像A上的投影p1在圖像B中的對應(yīng)點p2一定位于K2上；P在圖像B上的投影p2在圖像A中的對應(yīng)點p1一定位于K1上。由于幾何變換矩陣和其它計算的誤差，匹配點可能不會準確地出現(xiàn)在對應(yīng)的外極線上，而在外極線的鄰域內(nèi)。本發(fā)明方法在該鄰域內(nèi)搜索和外極線距離最近的點作為匹配點。由這幾組對應(yīng)點分別求出它們的三維坐標(biāo)。用直線段連接這些3D點，所得折線作為3D snake的初始位置。
snake是一條可變形的參數(shù)曲線c(s)＝(x(s)，y(s)，z(s))，s∈
。將該曲線離散化為N個點Ci(Xi，yi，Zi)(i＝1，2，…，N)，模型能量函數(shù)的離散表達式為 其中，Eint是內(nèi)部能量 Eint(i)＝α(d-|ci-ci-1|)+β|ci-1-2ci+ci+1|(7) 保證曲線變形過程中的連續(xù)性和光滑性。式中第一項約束使點平均分布，既能滿足連續(xù)性的要求，又不會產(chǎn)生曲線收縮的現(xiàn)象。在每次迭代結(jié)束時，點間的平均距離d的值被更新；第二項是曲率，使曲線保持平滑。α和β是權(quán)值。
外部能量Eext是保證模型收斂的外部力，吸引表示血管軸線的空間曲線移動到這樣的位置3D曲線在左右兩個圖像平面上的投影位于圖像上灰度最小、且灰度梯度最小的區(qū)域，即血管的2D中心線。根據(jù)式(1)和(5)，空間點P在圖像A和B上的投影坐標(biāo)(u1，v1)和(u2，v2)都可用該點的三維坐標(biāo)c＝(x1，y1，z1)來表示 [u1v1]T＝TL(c)，[u2v2]T＝TR(c) (8) 其中TL和TR都是以c為自變量的函數(shù)。Eext包含兩部分，分別對應(yīng)于兩個角度的二維投影 式中IL(TL(ci))＝IL(u1i，v1i)和IR(TR(ci))＝IR(u2i，v2i)分別是左右圖像點的灰度值； 和分別是左右圖像點的灰度梯度。由于造影圖像中，血管的灰度值比背景小，所以這里權(quán)重參數(shù)γ取正值。而且感興趣的是血管中心線，不是邊緣，所以λ也取正值。
(4)后續(xù)時刻血管軸線的四維重建 得到第一時刻血管段的三維軸線之后，對于其后的序列圖像，根據(jù)血管形狀不突變的特點，把前一時刻的三維血管軸線作為當(dāng)前時刻snake初始位置，實現(xiàn)對整個序列中血管軸線的四維重建。
在這一步驟中，本發(fā)明方法在snake能量函數(shù)中嵌入相似度匹配，通過對相鄰時刻之間的目標(biāo)進行配準，實現(xiàn)對血管的跟蹤。假設(shè)在足夠短的時間間隔內(nèi)，運動前后圖像中血管的灰度不發(fā)生變化，那么得到第一時刻的血管軸線之后，在開始對后續(xù)圖像進行運動跟蹤時，外部能量函數(shù)采用如下的形式 式中ILt和IRt分別表示時刻t左右圖像點的灰度值；ILt-1和IRt-1分別表示時刻t-1左右圖像點的灰度值；Δw是鄰域窗口的大小。(10)式中的前兩項與(9)式相同，表示吸引曲線向著圖像中灰度最小、灰度梯度最小的區(qū)域移動的圖像力。第三項是使得運動前后(時刻t-1和t)血管上對應(yīng)點的Δw×Δw鄰域內(nèi)的灰度變化最小的外部約束力。
對于求解snake能量函數(shù)最小化的問題，在計算機視覺發(fā)展的最初階段，一般都是采用變分法，迭代求解E-L偏微分方程。該方法不僅計算復(fù)雜，運算時間長，而且很容易收斂于局部最小值。采用動態(tài)規(guī)劃(“Using dynamic programming for solving variational problems invision，”IEEE Transactions on pattern analysis and machine intelligence.vol.12，no.9，pp.855-867，1990)可避免上述問題，由于迭代向著總能量減少的方向進行，當(dāng)總能量不再變化時，則迭代終止，因此保證了結(jié)果的收斂性。但是該方法的缺點是計算速度慢，需要較大的存貯空間。針對上述問題，本發(fā)明采用Williams(“A fast algorithm for active contours andcurvature estimation，”Computer Vision，Graphics and Image Processing，vol.55，no.1，pp.14-26，1992)提出的貪婪算法完成對最優(yōu)snake的搜索，該算法具有動態(tài)規(guī)劃方法的所有優(yōu)點，并且計算速度快，所需存貯空間小。
對于重建出的各時刻的三維血管軸線，根據(jù)式(8)即可算出其在左右兩個成像平面上的二維投影，即二維血管軸線。因此本發(fā)明方法將二維提取、三維重建和運動跟蹤集成到一個框架中完成，提高了運算精度和速度。
附圖4和5是對PHILIPS Integris CV單面全數(shù)字血管造影機臨床拍攝到的左冠造影圖像序列的實驗結(jié)果。采集速率為15幀/秒、圖像大小為512×512(像素)、灰階為256、像素尺寸為0.3mm。圖像序列的拍攝角度分別為RAO30°CAUD24°和LAO46°CRAN21°。圖4是第一時刻兩個角度的圖像對，在RAO圖像的左前降支(LAD)上選擇八個點，根據(jù)外極約束可找到其在LAO圖像中的匹配點。計算出它們的三維坐標(biāo)，用直線段連接這些3D點，得到用折線表示的3D snake的初始形狀，如圖4(a)所示。通過使能量函數(shù)((6)，(7)和(9)式)最小，snake不斷變形，最后得到第一時刻的3D血管軸線。采用三次B樣條曲線對這些離散3D點進行擬合，得到一條連續(xù)曲線表示的血管軸線，它在兩個成像平面上的投影就是左右圖像中的二維血管軸線，見圖4(b)(在原始圖像中用白色曲線表示)。
對于其后的序列圖像，把t時刻中模型的停留位置作為t+1時刻模型的初始位置，通過使能量函數(shù)((6)，(7)和(10)式)最小，模型在內(nèi)外力作用下變形，實現(xiàn)血管軸線序列的三維重建，即四維重建。圖5是選取的第4、6、8、10時刻的跟蹤結(jié)果，同樣在原始圖像上用白色曲線表示二維血管軸線。由此結(jié)果可見，當(dāng)心動周期結(jié)束時，即第十時刻的血管形狀與第一時刻的形狀(圖4(b))很相似，這也驗證了“在心動周期結(jié)束時，冠狀動脈血管會恢復(fù)到其在周期開始時的狀態(tài)”的結(jié)論。
權(quán)利要求
1、一種冠狀動脈血管軸線的四維重建方法，其特征是，首先采集兩個角度的覆蓋一個或多個心動周期的單面X射線冠狀動脈造影圖像序列、建立X射線血管造影系統(tǒng)在兩個角度的投影模型，并推導(dǎo)出兩幅不同角度的造影圖像之間的幾何變換關(guān)系，再對手動選取的第一時刻圖像中感興趣血管段的采樣點進行三維重建，并以連接各重建點所得折線作為snake模型的初始位置，通過snake變形獲得該段血管第一時刻的三維軸線，而對于圖像序列中的每一后續(xù)時刻，則均以其前一時刻的3D血管軸線作為snake模型的初始位置，通過snake變形獲得三維血管軸線，從而完成整個序列血管軸線的四維重建，具體步驟如下
a、以心電信號為同步信號，采集兩個角度的、覆蓋一個或多個心動周期的單面X射線冠狀動脈造影同步圖像序列，并對原始圖像進行初步的濾波、去噪、畸變校正處理，增強視覺效果；
b、建立X射線血管造影系統(tǒng)在兩個角度的投影模型，推導(dǎo)兩幅不同角度的造影圖像之間的幾何變換關(guān)系
設(shè)點s1和s2為兩次造影過程中X射線源焦點的位置，分別以s1和s2為原點，建立空間坐標(biāo)系s1x1y1z1和s2x2y2z2；坐標(biāo)系U1V1O1和U2V2O2為成像平面坐標(biāo)系；D1和D2分別是s1和s2到各自成像平面的垂直距離，空間血管上的點P在成像平面A和B上的投影點分別為p1(u1，v1)和p2(u2，v2)，
則s1x1y1z1和s2x2y2z2之間的幾何變換關(guān)系為
其中，R3×3為旋轉(zhuǎn)矩陣R＝RY(α2)·RX(β2)·RX(-β1)·RY(α1)
為平移向量
L1和L2為X射線源焦點s1和s2到旋轉(zhuǎn)中心的距離；(α1，β1)和(α2，β2)分別是A和B的造影角度；(x1，y1，z1)和(x2，y2，z2)分別表示空間點P在坐標(biāo)系s1x1y1z1和s2x2y2z2中的坐標(biāo)；
根據(jù)透視投影成像的幾何關(guān)系，空間點的三維坐標(biāo)與其二維投影坐標(biāo)之間的關(guān)系為
其中，
c、在序列中第一時刻的圖像中，通過人工取點獲得感興趣血管段在左右投影中的近似中心線，用折線表示，并利用外極約束得到兩個角度間對應(yīng)點的匹配；
d、對手動選取的采樣點進行三維重建；
e、連接各3D點，所得折線作為3D snake的初始位置，通過使預(yù)先設(shè)定的能量函數(shù)最小，snake在空間中變形，最終得到具有最小能量的最優(yōu)位置，就是第一時刻的3D血管軸線；
f、后續(xù)時刻血管軸線的四維重建
對于序列中的后續(xù)時刻，以前一時刻的3D血管軸線作為當(dāng)前時刻snake的初始位置，通過snake變形，完成整個序列中血管軸線的四維重建。
2、根據(jù)權(quán)利要求1所述的冠狀動脈血管軸線的四維重建方法，其特征是，所述X射線冠狀動脈造影圖像序列的兩個采集角度之間夾角的取值范圍為60°～120°。
3、根據(jù)權(quán)利要求1或2所述一種冠狀動脈血管軸線的四維重建方法，其特征是，在第一時刻圖像的感興趣血管段中選取的三維重建采樣點包括起點、終點和3～6個中間點。
全文摘要
一種冠狀動脈血管軸線的四維重建方法，屬醫(yī)學(xué)檢測技術(shù)領(lǐng)域。其技術(shù)方案是首先采集兩個角度的X射線冠狀動脈造影圖像序列、建立X射線血管造影系統(tǒng)在兩個角度的投影模型，并推導(dǎo)出兩個角度圖像之間的幾何變換關(guān)系，再對手動選取的第一時刻圖像中感興趣血管段的采樣點進行三維重建，并以連接各重建點所得折線作為初始位置，通過snake變形獲得該段血管第一時刻的三維軸線，而對于圖像序列中的每一后續(xù)時刻，則均以其前一時刻的3D血管軸線作為初始位置，通過snake變形獲得三維血管軸線，從而完成整個序列血管軸線的四維重建。本發(fā)明大大減少了操作者參與所引入的不確定性和誤差，提高了結(jié)果的可重復(fù)性，且操作簡便、效率高。
文檔編號G06T15/00GK101283911SQ200810055038
公開日2008年10月15日 申請日期2008年6月5日 優(yōu)先權(quán)日2008年6月5日
發(fā)明者正 孫 申請人:華北電力大學(xué)