孫 正

（華北電力大學電子與通信工程系, 河北 保定 071003）

目前臨床廣泛采用的診斷和治療冠心病的影像手段是 X 射線冠狀動脈造影（coronary angiography，CAG）和血管內超聲（intravascular ultrasound，IVUS）。CAG 通過靜止或動態(tài)觀察造影劑的充盈和消失情況來判斷血管解剖學形態(tài)異常的部位、性質和程度。該技術不僅能診斷缺血性心臟病及冠狀動脈畸形等疾病，而且對冠狀動脈內溶栓、PTCA（經皮腔內冠狀動脈成形術）等介入手術治療具有重要意義。但是，CAG只能反映血管腔被造影劑充填后的投影輪廓，因而不僅會存在盲區(qū)，而且不能提供血管壁的結構信息和病變程度。IVUS 是近年來應用于臨床診斷血管病變的一項新技術，將一微型高頻超聲探頭置于心導管頂端，利用心導管將超聲探頭導入血管腔內進行探測，再經電子成像系統(tǒng)來顯示血管組織結構和幾何形態(tài)的微細解剖信息?？筛鶕?jù)斑塊的聲學特征進行組織學分型，發(fā)現(xiàn)CAG 不能顯示的血管病變，觀察分叉處或血管重疊處的模糊病變。但是，由于采用16~40MHz 高頻探頭，影響了探測深度，只能對某一段病變血管進行測量，不能進入嚴重狹窄的管腔，并且無法確定截面的軸向位置和空間方向。

由于二維圖像的種種局限性，醫(yī)學圖像的三維重建和可視化已成為目前的研究熱點。基于多角度CAG 圖像的冠狀動脈三維重建[1]不僅能為醫(yī)生提供形象、直觀的三維血管圖像，而且可以對血管的有關參數(shù)（如長度、曲率等）進行定量測量。但是，該重建結果僅描述了血管內腔的空間位置和形態(tài)，無法得到管腔截面的形態(tài)和硬化斑塊的彌漫程度及組成。且重建時一般假設管腔橫截面為橢圓，事實上當發(fā)生狹窄時管腔形狀復雜多樣，狹窄多呈偏心型和不規(guī)則型，因此這種假設是不準確的。與CAG 相比，IVUS 有其獨特的圖像方位，可以顯示管腔的橫斷面、血管壁的形態(tài)結構、斑塊的組織形態(tài)學特征等，因此IVUS圖像的三維重建有利于評價沿血管長軸方向管腔的邊界、截面輪廓以及斑塊的幾何形態(tài)等。傳統(tǒng)的IVUS 三維重建是把一系列的超聲圖像按采集順序疊加起來形成一個三維直血管段[2]。由于IVUS 本身不能提供每一幀圖像的空間幾何信息，因此這種方法沒有考慮在圖像獲取過程中導管的彎曲和扭曲，其結果也是不準確的。

從上面的分析可知，基于CAG 圖像的血管三維重建和傳統(tǒng)的IVUS 三維重建各有優(yōu)缺點，而且具有互補性。本文通過將由IVUS 圖像獲取的血管橫截面信息，和由基于造影圖像的三維重建獲得的導管空間幾何信息結合起來，克服IVUS 和CAG 各自顯示血管形態(tài)的局限性，準確重建出血管的解剖結構和反映血管的真實彎曲和扭曲，得到病變的準確位置和形態(tài)，輔助冠心病的診治。

1 融合方法

基于CAG 與IVUS 融合的血管3D 重建主要包括4 個步驟：① 圖像采集；② 從CAG 圖像中提取超聲導管的回撤路徑，并進行三維重建；③ 對IVUS 圖像序列進行分割，提取血管壁內外膜的邊緣；④ 將基于CAG 三維重建獲得的導管空間幾何信息和由IVUS 獲得的血管橫截面信息進行數(shù)據(jù)融合，完成重建。

1.1 圖像采集

在IVUS 的操作過程中，造影和IVUS 分別同步顯示探頭在管腔內的部位和相應血管壁的結構形態(tài)（如圖1 所示，CAG 中的箭頭指明獲得IVUS 圖像的位置）。在CAG 中發(fā)現(xiàn)可疑血管后，向目標血管內注入硝酸甘油，在X 射線透視圖像的指導下將超聲導管穿越病變部位，到達血管遠端，將超聲探頭與超聲成像儀連接去除偽影后，經馬達控制勻速回撤導管，并記錄圖像。由于心臟的周期性運動和呼吸的影響，很難獲得對應于同一時刻的圖像序列，因此本文采用ECG（心電）門控的方法，在相同心臟相位處采集圖像。

對于CAG 圖像的采集，本文采用僅在導管回撤路徑的起點拍攝一對造影圖像的方法，從而減少了X 射線的劑量，提高了臨床應用價值。實驗中采用機械式超聲探頭，超聲換能器位于一可彎曲的軸心頭端，軸心在外鞘管內以1800 轉/分的速度旋轉，而鞘管是固定不動的，因此可以保證回撤路徑的穩(wěn)定，即在回撤導管的過程中，帶有超聲換能器的導管尖端的回撤軌跡不會偏離在回撤路徑起點拍攝的造影圖像中顯示的導管影像。成像過程中記錄造影角度和X 射線源焦點至接收屏的距離。

1.2 造影圖像中導管路徑和管腔邊緣的提取和三維重建

造影圖像中導管路徑的提取和重建類似于血管軸線的提取和重建。有些算法是直接把血管軸線作為導管路徑的近似[3]，事實上二者并不重合，如圖2 所示。文獻[4]在建立了CAG 系統(tǒng)在兩個近似垂直角度的透視投影成像模型（圖3）之后，根據(jù)在造影過程中同步記錄的距離和角度參數(shù)，得到成像系統(tǒng)的幾何變換矩陣。然后利用3D snake 模型，表示血管軸線的曲線直接在空間中變形，完成血管軸線的三維重建。該方法將圖像中血管的二維軸線提取和三維重建集成到一個框架中完成，避免了基于外極約束的兩個角度間的逐點匹配，提高了重建精度和運算速度。因此，本文采用類似方法完成導管路徑的三維重建，3D snake 的初始位置采用手動取點獲得，即在導管路徑的一個投影上手動選取若干采樣點（一般選取回撤路徑的起點、終點和3～4 個中間點即可），然后根據(jù)外極約束得到這些點在另一投影上的對應點。由這幾組對應點分別求出它們的三維坐標[1]，用直線段連接這些3D 點，所得折線作為3D snake 的初始位置。

圖2 導管路徑與血管中心線示意圖

按照成像系統(tǒng)的幾何變換矩陣，將3D 導管路徑向兩個成像平面反投影，即可得到對應的2D路徑。對于2D 路徑上的每個點，通過在垂直于路徑的方向上，尋找灰度梯度的兩個極大值，即可完成對血管腔左右邊緣的提取（圖4）。之后，在假設管腔的橫截面是橢圓的前提下，可完成整個管腔的三維重建[5]。

圖3 造影系統(tǒng)在兩個角度的成像示意圖

圖4 造影圖像中管腔邊緣的提取

1.3 血管內超聲圖像的邊緣提取

正常血管壁在超聲圖像上呈現(xiàn)回聲不同的三層環(huán)形結構[6]：內層為一強回聲亮環(huán)，在組織學上由內膜及內彈性膜組成；中層相當于中膜，呈現(xiàn)為低回聲暗區(qū)；外層為較亮的強回聲帶，由外彈性膜和外膜構成。序列中前后幀之間非常相似，具有很強的相關性。

目前臨床分割IVUS 圖像常采用手動方法，由經驗豐富的醫(yī)生用鼠標驅動畫筆逐一在每幀中繪出斑塊―中膜、中膜―外膜以及斑塊―內腔的邊界線。該方法不僅對操作者的專業(yè)知識和技能要求較高，而且結果的可重復性差，主觀性較強。近年來，提出了很多計算機輔助的分割方法，為了保證結果的精度，仍然需要操作者一定程度的手動參與。大量研究表明，對于IVUS 圖像，結合動態(tài)規(guī)劃等處理方法的snake 模型能夠獲得良好的邊緣提取能力[7]。操作者只需在首幀中手動選擇目標輪廓上的幾個特征點，連接這些點所形成的多邊形作為snake 的初始位置。對于后續(xù)幀，將前一幀的提取結果作為下一幀snake 的初始位置，可完成對連續(xù)多幀IVUS 圖像的分割，大大節(jié)省計算時間。對一幀IVUS 圖像的分割結果如圖5 所示。

圖5 完成輪廓提取的IVUS 圖像

1.4 IVUS 與CAG 的融合

在獲得血管的2D 橫截面信息和導管回撤路徑的3D 空間幾何信息之后，就需要將在各個心臟相位上采集到的超聲圖像沿回撤路徑準確排列起來。主要需解決兩個問題：確定各IVUS 幀的3D 軸向位置和空間方位（包括相鄰幀間的相對方位和各幀的絕對方位）。

1.4.1 超聲圖像三維軸向位置的確定

在對感興趣血管段進行IVUS 檢查的過程中，將超聲導管插入管腔內后，推送至遠端，然后經馬達控制勻速地從遠端向近端連續(xù)拉出導管，進行超聲檢測，順序等距采集圖像?？梢哉{節(jié)回撤速度，從而根據(jù)需要調節(jié)切面間距。采用CAG 圖像重建出導管軸線之后，根據(jù)已知的切面間距依軸向將各幀IVUS 圖像順序排列，即可確定出各幀的軸向位置。

1.4.2 相鄰幀超聲圖像間相對方位的確定

確定各相鄰IVUS 幀間的相對方位，也就是考慮導管在回撤過程中的扭轉問題。當導管的3D回撤路徑不是平面曲線時，如圖6 所示，是沿著一個立方體邊沿的回撤軌跡，A 和B 點分別表示回撤路徑的起點和終點，路徑上的箭頭表示在該點采集到的超聲圖像的方向，也即導管在該點所在平面的法線方向。顯然當導管從一個平面進入另一個平面時，超聲圖像的方向就會發(fā)生變化。而且，當導管分別沿實線和虛線表示的兩條路徑回撤時，盡管起點和終點相同，所采集的超聲圖像的方向也是不同的：對于實線軌跡，在A 和B處的方向正好相差180°；而虛線軌跡的起點和終點處的方向是相同的。如果在重建過程中沒有考慮此問題，則可能造成實際的斑塊在血管彎曲處的外側，而重建出的斑塊位于彎曲處內側的現(xiàn)象，這兩種情況下的血液動力學特性是完全不同的。

圖6 導管回撤過程中的扭曲示意圖

本文采用建立微分幾何模型的方法解決該問題，即在重建后的3D 導管路徑上建立各幀超 聲圖像的局部坐標系。設c(s)(0 1)s≤ ≤ 是完成 B 樣條曲線擬合后的3D 導管路徑，則在c(s)上某點c(s0)處的局部正交坐標系如圖7 所示，該坐標系是Frenet-Serret 標架，其原點為c(s0)，3 個坐標軸分別為單位切矢量t、單位主法矢量n 和單位副法矢量b，它們可由曲線方程得到[8]。對于沿導管路徑軸向分布的各幀圖像均建立Frenet-Serret 標架，即可確定各幀間的相對方位，如圖8 所示。

圖7 導管路徑上各幀IVUS 圖像的局部坐標系

圖8 超聲圖像序列中相鄰幀間相對方位的確定

1.4.3 超聲圖像序列絕對方位的確定

在確定了相鄰幀間的相對方位后，各幀的絕對方位仍是未知的。這就好比在腿上穿襪子，當腿保持不動時，襪子能夠圍繞腿自由旋轉，但是僅存在一個方位是最適合的。這里導管路徑就如同腿，而IVUS 序列則如同襪子，因此確定絕對方位就如同尋找襪子最適合的方位。

首先利用圖8 所示的局部坐標系確定各幀的初始方位。在IVUS 圖像中，導管位于圖像中心，分割出的血管壁輪廓的重心一般不與導管重合，如圖9 所示，其中C 點表示導管，OC為橢圓輪廓的中心（即在假設血管橫截面是橢圓時，基于CAG 的三維重建中所對應的血管中心線位置），OI為從超聲圖像中提取出的血管壁輪廓的重心

由于血管中心線和導管路徑不重合，在血管同一位置處的超聲圖像輪廓和基于CAG 重建出的橢圓輪廓方位不一致，把橢圓輪廓投影到對應的超聲圖像上，如圖9 所示。同樣定義橢圓輪廓的離心向量μ 來表示血管軸線偏離導管的程度

超聲圖像的匹配誤差可用向量ρ 的模ε 和ρ 與μ的夾角θ 表示，如圖9 所示。第i 幀超聲圖像的εi和θi為

圖9 超聲圖像偏心距離和偏心夾角示意圖

利用這些數(shù)值，在每一個窗口位置處，計算可靠性權重因子

并將其應用于各IVUS 幀，從而獲得各幀圖像的絕對方位。

1.5 重建結果的顯示

在對IVUS 圖像序列完成邊緣提取并確定各幀的空間位置后，利用基于NURBS 表面擬合的表面提取法[5]完成血管的重建，即管腔內外表面的繪制。此后，可利用虛擬現(xiàn)實造型語言（Virtual Reality Modeling Language，VRML）來顯示重建結果。

2 實驗結果及分析

2.1 實驗結果

為了驗證算法的可行性，本文采用臨床采集 的CAG 和IVUS 數(shù)據(jù)進行了實驗。編程所用的計算機配置為Intel P4 2.8GHz CPU、1GB 內存。編程環(huán)境為Visual C++6.0，其中部分數(shù)學計算采用Matlab 6.5 完成。

CAG 圖像是在型號為Philips Integris CV 的全數(shù)字單面X 射線血管造影機上采集到的。采用Jomed Endosonic 超聲成像儀，探頭為2.9F Jomed單軌機械探頭，頻率為30MHz。在回撤導管的過程中，探頭導管以1800 轉/分作360°旋轉時可連續(xù)獲得30 幀/秒的血管橫軸實時切面圖像。圖10是右冠CAG 和IVUS 圖像采集示意圖。圖11 是在導管回撤路徑的起點拍攝的一對造影圖像，拍攝角度分別為RAO75o和LAO20o，圖中提取出了血管腔的左右輪廓，用白色線條表示，黑色線條表示導管路徑。圖12 是任選的一幀IVUS 圖像，其中提取出了血管壁內外膜的輪廓。圖13 是將CAG和IVUS數(shù)據(jù)融合后得到的血管段三維重建結果。

圖10 右冠CAG 和IVUS 圖像采集示意圖

圖11 在導管回撤路徑起點拍攝的一對造影圖像

圖12 IVUS 圖像

圖13 重建結果

2.2 誤差分析

由于很難得到真實人體血管管腔的形態(tài)，因此對重建結果精度的定量分析也很困難。本文在后續(xù)工作中考慮采用模型的CAG 和IVUS 圖像數(shù)據(jù)進行實驗，定量估計算法的誤差。本節(jié)僅對重建過程中可能產生的誤差進行定性分析。從重建步驟可知，結果的精度取決于以下幾方面：

（1） 圖像采集的質量

· 采集時間：單面CAG 系統(tǒng)在兩個角度拍攝的圖像是在不同時刻獲得的，由同步記錄的ECG 信號可以控制造影時間，使得兩幅圖像盡量同步。同時由于心臟的劇烈運動，超聲圖像的采集也需要采用ECG 門控的方法。但ECG 信號本身也存在誤差，而且病人的呼吸和身體的移動也會使血管的空間位置發(fā)生變化。在實際操作過程中可通過讓病人屏住呼吸，提高采集速度，減小該影響。

· 造影角度：造影角度的準確性直接影響重建精度。臨床應用中，通常由操作人員根據(jù)心臟的形態(tài)、病人狀態(tài)和醫(yī)生的觀察習慣等選擇角度。 由于機械裝置的限制， 一般在LAO90o~RAO90o和CAUD45o~ CRAN45o的范圍內選擇。同時為了進行3D 重建，要求左右角度之差應在60o~120o范圍內[9]。

（2） IVUS 圖像分割和提取的質量

目前超聲圖像分割的方法主要有全手動、半自動和全自動3 種。全手動方法可靠性較高，但缺點是對操作人員的技術水平和專業(yè)知識要求較高，結果不可避免的受到操作者技術和主觀因素的影響，可重復性差，工作量大，對處理大量數(shù)據(jù)不現(xiàn)實?，F(xiàn)有全自動算法僅適用于圖像質量較高的情形，當圖像受污染較嚴重，或由于病變造成管腔截面形狀比較復雜時，很難得到滿意結果。半自動的方法包括兩種：自動獲取近似形狀，再手動修正；手動粗略勾畫，再自動修正。前者采用簡單的圖像處理技術獲取初始輪廓，手動修正時，操作者不僅要利用圖像本身的特征信息，而且還要結合局部解剖和病理知識，同樣很耗時，可重復性差。后者減少了操作者的參與，僅需手動勾畫輪廓的粗糙形狀，后續(xù)自動提取過程會對其進行修正，因此操作者的參與對最終結果的影響是間接的。在參數(shù)設置相同的情形下，初始輪廓一定程度上的變化不會影響最終結果，因此可重復性很高。本文采用的snake 模型方法即屬此類。同時與其它非模型方法相比，它將輪廓看作一條連續(xù)曲線，即使在相應圖像特征很弱甚至缺失的情形下，仍可提取出連續(xù)的輪廓。

（3） 導管回撤路徑三維重建的精度

首先，重建精度受到CAG 成像系統(tǒng)幾何參數(shù)的影響。由于機械制造的限制，造影系統(tǒng)提供的角度和距離參數(shù)存在一定誤差。由于成像系統(tǒng)的標定過于復雜，需要特殊的實驗裝置，所以臨床一般不進行標定。利用血管樹分叉點優(yōu)化幾何變換矩陣的方法，可在一定程度上降低系統(tǒng)參數(shù)誤差的影響[9]。此外有些情況下，在CAG 圖像中并不能清晰觀察到導管的影像，甚至完全不可見。此時只能采用血管腔的中軸線作為回撤路徑的近似。但是，如圖2 所示，當導管在回撤過程中偏離管腔軸線時，重建出的3D 血管軸線的總長通常會大于真實的導管路徑總長。

（4） IVUS 幀的定位精度

由于IVUS 本身不能提供各幀的空間信息，因此將各幀圖像映射到導管路徑的正確位置和方向上是準確重建血管的關鍵。其中軸向位置的確定較為簡單，只要保證勻速回撤和準確記錄回撤起點即可?？臻g方位的確定比較困難，首先要考慮導管的扭轉問題，即確定各幀的相對方位。其次需確定各幀的絕對方位，為每幀找到其最合適的位置。目前國外各研究組均致力于尋找解決此問題的最佳方法，這也是作者今后的主要研究方向。

2.3 運算時間和自動化程度

為了提高融合系統(tǒng)的臨床實用性，所采用的方法應遵循兩個原則：一是盡量減少對臨床手術的干擾；二是盡量提高重建精度。如前所述，本文算法的基本步驟包括圖像采集、圖像分割、重建導管路徑和融合。根據(jù)上述原則，只有第一步在心導管室完成，即在進行介入檢查或治療中完成圖像采集和相關參數(shù)的記錄。后續(xù)步驟均可在手術之外的其它時間進行，不必要求同步實時處理。其中導管路徑的三維重建和IVUS 序列的分割均采用基于snake 模型的半自動方法，僅在獲取snake 初始位置時（對于超聲序列，僅在獲得首幀初始snake 時）需手動參與，后續(xù)步驟均自動完成。對于snake 能量函數(shù)最優(yōu)值的求解，本文采用Williams 快速算法[10]，較變分法和動態(tài)規(guī)劃法而言，可大大縮短計算時間[11]。在本文實驗采用的計算機配置下，對于長度為90mm 的回撤路徑，重建路徑僅需約0.2 秒；分割20 幀超聲圖像，需要約1.1 秒；完成各幀超聲圖像的空間定位需要約1.4 秒。若采用更高配置的計算機，計算時間將會進一步縮短。

3 結 論

X射線血管造影和血管內超聲在顯示血管形態(tài)上存在不同的局限性，但是由于這兩種技術具有互補性，因此通過將它們融合，可以克服各自的不足，得到準確的血管空間位置和形態(tài)。本文提出了基于IVUS 和CAG 數(shù)據(jù)融合的三維重建血管的方法。在完成兩種圖像的采集后，采用結合動態(tài)規(guī)劃的snake 模型對各幀IVUS 圖像進行分割，提取出血管壁內外膜的輪廓。然后采用基于3D snake 的半自動方法從CAG 圖像中重建超聲導管的回撤路徑。在確定各幀IVUS 圖像在導管路徑軸向上的位置后，通過Frenet-Serret 標架確定各相鄰IVUS 幀之間的相對方位，利用非迭代的統(tǒng)計優(yōu)化方法確定各幀超聲圖像的絕對方位，完成IVUS 和CAG 的融合。采用臨床數(shù)據(jù)驗證了方法的可行性。

進一步的研究是改進確定IVUS 幀空間方位的算法，提高融合精度。考慮采用模型進行實驗，定量估計算法的精度。同時考慮對重建結果進行四維顯示，通過將三維融合結果和血管的運動相結合，準確地顯示血管的形態(tài)、方位、運動情況和血管橫截面信息等，為冠心病的診治提供依據(jù)。

[1] 郁道銀, 黃家祥, 陳曉冬, 等. 基于B 樣條的冠狀動脈樹骨架三維重建方法[J]. 工程圖學學報, 2005, 26(2): 95-100.

[2] Birgelen C, Vrey EA, Mintz GS, et al. ECG-gated three-dimensional intravascular ultrasound: feasibility and reproducibility of the automated analysis of coronary lumen and atherosclerotic plaque dimensions in humans [J]. Circulation, 1997, 96(9): 2944-2952.

[3] Prause GPM, DeJong SC, McKay CR, et al. Towards a geometrically correct 3D reconstruction of tortuous coronary arteries based on biplane angiography and intravascular ultrasound [J]. International Journal of Cardiac Imaging, 1997, 13(6): 451-462.

[4] 孫 正, 郁道銀, 姜 浩. 基于變形模型的心血管三維運動跟蹤[J]. 光電工程, 2006, 33(6): 24-27.

[5] 陳曉冬, 黃家祥, 郁道銀, 等. 基于造影圖像的冠狀動脈血管表面三維重建方法[J]. 工程圖學學報, 2005, 26(3): 111-116.

[6] Cardinal MR, Meunier J, Soulez G, et al. Intravascular ultrasound image segmentation: a three-dimensional fast-marching method based on gray level distributions [J]. IEEE Transactions on Medical Imaging, 2006, 25(5): 590-601.

[7] Luo Z, Wang Y, Wang W. Estimating coronary artery lumen area with optimization-based contour detection [J]. IEEE Transactions on Medical Imaging, 2003, 22(4): 564-566.

[8] 朱心雄. 自由曲線曲面造型技術[M]. 北京: 科學出版社, 2000. 1-14.

[9] 黃家祥, 郁道銀, 陳曉冬, 等. 冠脈樹三維重建中幾何變換矩陣的優(yōu)化[J]. 中國生物醫(yī)學工程學報, 2005, 24(2): 189-193.

[10] Williams DJ, Shah M. A fast algorithm for active contours and curvature estimation [J]. Computer Vision, Graphics and Image Processing, 1992, 55(1): 14.

[11] 孫 正, 郁道銀, 姜 浩. 造影圖像序列中血管運動跟蹤的Snake 方法[J]. 中國醫(yī)學影像技術, 2004, 20(11): 1779-1783.