定量分析左心室射血運(yùn)動(dòng)特征的追蹤方法設(shè)計(jì)

鄭南南，黃 鋼，林劍圣

（1.上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海 200093；2.上海健康醫(yī)學(xué)院，上海 201318）

0 引言

心血管疾?。–ardiovascular Disease，CVD）是全球最常見的非傳染性疾?。?］。我國(guó)心血管疾病發(fā)生率與死亡率正處于上升階段，其中死亡率位居各類疾病首位［2］，如何有效預(yù)防和診療心血管疾病成為重大公共衛(wèi)生問題。

左心室射血運(yùn)動(dòng)為全身輸送營(yíng)養(yǎng)物質(zhì)，是心臟功能的重要組成部分，異常的左心室運(yùn)動(dòng)提示可能患有心血管疾病，如冠脈梗塞引起的心肌供血不足會(huì)導(dǎo)致左心室射血無力。先進(jìn)的醫(yī)學(xué)成像技術(shù)，如心臟超聲（Cardiac Ultra?sound，CU）、計(jì)算機(jī)斷層掃描（Computed Tomograph，CT）以及心臟電影磁共振成像（Cardiac Cine Magnetic Resonance Imaging，CCMRI）等為醫(yī)師提供了動(dòng)態(tài)觀察左心室是否正常運(yùn)動(dòng)的手段。然而，僅靠觀察心臟運(yùn)動(dòng)狀態(tài)診斷心血管疾病尚顯不足，科研工作者正致力于準(zhǔn)確量化心臟圖像中左心室運(yùn)動(dòng)的研究。Valizadeh 等［3］在形變配準(zhǔn)中引入心內(nèi)膜形狀信息以提高配準(zhǔn)精度，但該法僅測(cè)量了左心室徑向應(yīng)變，對(duì)左心室完整運(yùn)動(dòng)的量化不夠充分；Leong 等［4］結(jié)合釓增強(qiáng)和標(biāo)記的磁共振圖像融合信息分析心肌梗死患者的左心室運(yùn)動(dòng)，然而被標(biāo)記的磁共振圖像時(shí)空分辨率不高，較少應(yīng)用于臨床診斷心血管疾??；Mantilla 等［5］通過字典學(xué)習(xí)方法對(duì)左心室正常與異常運(yùn)動(dòng)進(jìn)行分類，但該法不能導(dǎo)出左心室的運(yùn)動(dòng)細(xì)節(jié)，如位移場(chǎng)信息等；Qiao 等［6］分別采用傳統(tǒng)配準(zhǔn)方法與新興卷積神經(jīng)網(wǎng)絡(luò)（Convolutional Neural Network，CNN）追蹤左心室運(yùn)動(dòng)，與傳統(tǒng)配準(zhǔn)方法相比，當(dāng)訓(xùn)練數(shù)據(jù)與測(cè)試數(shù)據(jù)來自相同健康志愿者時(shí)，CNN表現(xiàn)出更高的運(yùn)動(dòng)估計(jì)精度。然而CNN 評(píng)估心臟病患者左心室運(yùn)動(dòng)的準(zhǔn)確度并不高，需要加入更多患者數(shù)據(jù)加強(qiáng)其泛化能力。目前基于像素灰度學(xué)習(xí)的左心室運(yùn)動(dòng)追蹤比較困難，因?yàn)槿绻豢紤]模擬數(shù)據(jù)，真實(shí)的左心室運(yùn)動(dòng)場(chǎng)幾乎不可能獲取。Tuyisenge 等［7］利用非線性光流法恢復(fù)左心室二維運(yùn)動(dòng)場(chǎng)，由于左心室心肌具有大尺度形變的特點(diǎn)，光流估計(jì)中的小運(yùn)動(dòng)假設(shè)并不符合左心室的運(yùn)動(dòng)事實(shí)。

針對(duì)以上問題，本文設(shè)計(jì)了一個(gè)通用左心室運(yùn)動(dòng)追蹤框架，可用于不同心臟成像模態(tài)。該框架不需要對(duì)各類心臟疾病數(shù)據(jù)進(jìn)行學(xué)習(xí)，能依據(jù)左心室輪廓特征，快速恢復(fù)像素級(jí)左心室位移場(chǎng)，詳細(xì)展示左心室的運(yùn)動(dòng)細(xì)節(jié)。此外，該框架還可根據(jù)求得的位移場(chǎng)信息，分析左心室射血運(yùn)動(dòng)的各種量化參數(shù)，以便多角度輔助診斷心血管疾病。左心室運(yùn)動(dòng)追蹤初始工作建立在心臟電影磁共振成像基礎(chǔ)上，這是由于其具有無輻射、軟組織對(duì)比度高、成像角度豐富等獨(dú)特優(yōu)勢(shì)，廣泛應(yīng)用于臨床。考慮到左心室輪廓豐富的形狀屬性能夠?yàn)檫\(yùn)動(dòng)追蹤提供可靠特征［8］，本文首先利用三角網(wǎng)格構(gòu)建左心室表面模型，該模型是輪廓點(diǎn)連通關(guān)系的表達(dá)；然后擬合輪廓點(diǎn)局部曲面，計(jì)算輪廓點(diǎn)彎曲能量；其次，基于映射準(zhǔn)則估計(jì)所有當(dāng)前幀輪廓點(diǎn)最可能的下幀位置，恢復(fù)初始輪廓位移場(chǎng)；最后，平滑優(yōu)化輪廓位移場(chǎng)并計(jì)算左心室的位移和形變參數(shù)。

1 實(shí)驗(yàn)方法

擬建立一個(gè)基于心臟電影磁共振短軸圖像的左心室運(yùn)動(dòng)追蹤方法框架，恢復(fù)整個(gè)心動(dòng)周期的左心室輪廓位移場(chǎng)。如圖1 所示，整個(gè)方法流程主要包括4 個(gè)階段：數(shù)據(jù)預(yù)處理、左心室表面模型生成、特征追蹤和位移場(chǎng)優(yōu)化。

1.1 數(shù)據(jù)預(yù)處理

數(shù)據(jù)預(yù)處理包括左心肌分割、輪廓平滑和插值。左心肌分割的目的為提取心內(nèi)膜與心外膜輪廓，是本文的主要研究對(duì)象。目前常用的左心肌分割方法為圖像驅(qū)動(dòng)、主動(dòng)輪廓模型和CNN［9-11］。本文采用CVI 心臟分析軟件中左心肌分割功能獲取的二值輪廓圖像作為后續(xù)流程的基礎(chǔ)圖像。

Fig.1 Framework of motion tracking algorithm of left ventricle圖1 左心室運(yùn)動(dòng)追蹤方法框架

左心肌分割會(huì)因各類噪聲，如圖像質(zhì)量差、左心肌周圍組織干擾而得到不平滑輪廓。為確保運(yùn)動(dòng)追蹤結(jié)果的可靠性，應(yīng)消除局部尖銳的輪廓形狀。為此，利用2 階高斯平滑算法［12］，使輪廓在整個(gè)心動(dòng)周期的變化均趨于平滑。

由于覆蓋心臟的電影磁共振短軸圖像層間分辨率大于層內(nèi)分辨率，因此在平滑輪廓后需要對(duì)短軸圖像進(jìn)行層間插值，以構(gòu)造更為精細(xì)的左心室表面模型，有利于左心室的三維運(yùn)動(dòng)追蹤。本文采用基于輪廓形狀的插值算法［13-14］，以彌補(bǔ)層間分辨率的不足。

1.2 左心室表面模型生成

將左心室各層面輪廓按照心尖到心底的順序堆疊，并采用三角網(wǎng)格的方式連接相鄰輪廓，由此生成左心室表面模型。狄洛尼三角剖分具有空?qǐng)A性、最大化最小角等優(yōu)良特性［15］。基于狄洛尼三角剖分思想，分兩步構(gòu)造相鄰輪廓的三角網(wǎng)格：第一步是確定一組對(duì)稱最近鄰點(diǎn)對(duì)，該點(diǎn)對(duì)中任一輪廓點(diǎn)的歐式距離最鄰近點(diǎn)均為點(diǎn)對(duì)中的另一輪廓點(diǎn)，如圖2 左所示；第二步是以某個(gè)對(duì)稱點(diǎn)對(duì)作為起始三角網(wǎng)格的兩個(gè)頂點(diǎn)P1和P2，比較點(diǎn)P1與P4、點(diǎn)P2與P3的歐式距離，選擇歐式距離更小的輪廓點(diǎn)作為第3 個(gè)頂點(diǎn)，至此一個(gè)三角網(wǎng)格構(gòu)造完畢，如圖2 右所示。該方法能構(gòu)造出使所有三角網(wǎng)格邊長(zhǎng)總和最小的左心室表面模型。

Fig.2 Steps of construction of triangular mesh圖2 三角網(wǎng)格構(gòu)造步驟

對(duì)每個(gè)相鄰輪廓均實(shí)施上述構(gòu)網(wǎng)法則，便可得到左心室表面模型，如圖3 所示。該模型能表達(dá)輪廓的局部形狀屬性，包括輪廓點(diǎn)與其鄰域點(diǎn)的連通關(guān)系信息，為后續(xù)輪廓特征提取奠定了基礎(chǔ)。

Fig.3 Surface model of left ventricle original contour（a）and interpolated contour（b）圖3 左心室表面模型原始輪廓（a）與插值輪廓（b）

1.3 特征追蹤

輪廓點(diǎn)與其鄰域點(diǎn)形成的自然連通域是輪廓點(diǎn)局部形狀屬性的表征，如圖4（a）所示。然而，自然連通域只是由有限個(gè)輪廓點(diǎn)構(gòu)成的離散曲面。為了能較好地提取輪廓點(diǎn)局部形狀屬性，應(yīng)對(duì)自然連通域進(jìn)行局部曲面擬合。選擇二元二次函數(shù)作為曲面擬合函數(shù)［16］，表示為：

式中，x、y、z 為輪廓點(diǎn)坐標(biāo)，a1～a5為需要求解的5 個(gè)逼近系數(shù)。采用最小二乘估計(jì)，以矩陣的形式表示，求得使逼近誤差最小的系數(shù)向量a，表示為：

式中，A為n×5的矩陣，第i行每列元素分別為，n為構(gòu)成自然連通域的輪廓點(diǎn)數(shù)目，z為輪廓點(diǎn)坐標(biāo)值z(mì)的列向量。

計(jì)算出系數(shù)向量a后，便可以得到擬合的局部曲面，如圖4（b）所示。根據(jù)求得的曲面函數(shù)，分析其曲面微分性質(zhì)，選擇最大主曲率k1和最小主曲率k2作為形狀屬性參數(shù)［17］。

Fig.4 Local surface fitting discrete surface（a）and continuous surface（b）圖4 局部曲面擬合離散曲面（a）與連續(xù)曲面（b）

左心室運(yùn)動(dòng)會(huì)使輪廓產(chǎn)生彎曲變形，而主曲率能夠度量輪廓的彎曲變形程度。將輪廓看作薄彈性板，彈性板的形變勢(shì)能不受旋轉(zhuǎn)與平移運(yùn)動(dòng)影響，僅由主曲率和彈性板的材料常數(shù)決定。具體表示為：

式中，εpo為輪廓點(diǎn)局部曲面勢(shì)能，k1和k2分別為最大與最小主曲率，M為材料常數(shù)。

將輪廓點(diǎn)局部曲面勢(shì)能映射到其鄰域范圍，得到輪廓的勢(shì)能分布情況，如圖5 所示。由圖5 可知，左心室輪廓的整體勢(shì)能較低，舒張期輪廓?jiǎng)菽鼙仁湛s期更低，說明舒張期輪廓的形變程度更小。

Fig.5 Potential energy distribution of contours diastole period（a）and systolic period（b）圖5 輪廓?jiǎng)菽芊植际鎻埰冢╝）與收縮期（b）

將輪廓點(diǎn)的局部曲面勢(shì)能作為形狀特征，根據(jù)映射準(zhǔn)則進(jìn)行特征追蹤。該映射準(zhǔn)則以輪廓在很短時(shí)間間隔內(nèi)勢(shì)能變化微小的假設(shè)為前提，表示為：

式中，kp1和kp2分別為當(dāng)前幀輪廓點(diǎn)的最大與最小主曲率，kf1和kf2分別為下幀輪廓點(diǎn)的最大與最小主曲率，Pm為當(dāng)前幀輪廓點(diǎn)的下幀匹配點(diǎn)，該匹配點(diǎn)為當(dāng)前幀輪廓點(diǎn)位移范圍內(nèi)的最佳映射點(diǎn)。采用對(duì)稱最近鄰法確定當(dāng)前幀輪廓點(diǎn)的下幀匹配范圍［18］。

1.4 位移場(chǎng)優(yōu)化

將匹配準(zhǔn)則應(yīng)用于每個(gè)當(dāng)前幀輪廓點(diǎn)，得到輪廓初始位移場(chǎng)，如圖6（a）所示。然而，單憑孤立的點(diǎn)到點(diǎn)匹配過程恢復(fù)的初始位移場(chǎng)在空間上沒有局部相干性和整體協(xié)調(diào)性，這是由于當(dāng)前幀輪廓點(diǎn)的匹配不受其他鄰域輪廓點(diǎn)的約束，是一個(gè)空間獨(dú)立過程，獨(dú)立匹配得到的初始位移場(chǎng)并不符合左心室平滑的運(yùn)動(dòng)特征。因此，為了得到合理準(zhǔn)確的位移場(chǎng)，需要對(duì)初始位移場(chǎng)進(jìn)行平滑優(yōu)化，表示為：

式中，μ和λ為平滑因子，d(p)ini為當(dāng)前幀輪廓點(diǎn)初始位移向量，d(p)ave為鄰域輪廓點(diǎn)的平均位移向量，d(p)smo1和d(p)smo2為平滑優(yōu)化的位移向量。對(duì)每個(gè)當(dāng)前幀輪廓點(diǎn)進(jìn)行一定次數(shù)的兩輪平滑迭代，便可以得到優(yōu)化后的位移場(chǎng)，如圖6（b）所示。

Fig.6 Displacement field of left ventricle original displacement field（a）and optimal displacement field（b）圖6 左心室位移場(chǎng)原始位移場(chǎng)（a）與優(yōu)化位移場(chǎng)（b）

2 實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)對(duì)象包括35 例女性和66 例男性，年齡14-88 歲，心率46～125 次/min，其心臟電影磁共振短軸圖像數(shù)據(jù)均由GE1.5T 掃描儀通過穩(wěn)態(tài)自由進(jìn)動(dòng)（Steady State Free Preces?sion，SSFP）序列獲取。采集參數(shù)：圖像矩陣大小為256×256，視野為360mm×360mm，重復(fù)時(shí)間為3.5ms，回波時(shí)間為1.4ms，層厚為6～7mm，層間距為7～10mm，可用于左心室運(yùn)動(dòng)追蹤的幀數(shù)為20～28，層數(shù)為3～9。

2.2 實(shí)驗(yàn)結(jié)果

根據(jù)心臟電影磁共振短軸圖像數(shù)據(jù)，首先計(jì)算101 例受試者的左心室射血分?jǐn)?shù)（Ejection Fraction，EF），將其分為射血無力組（46 例，EF 小于50%）和射血正常組（55 例，EF大于50%），其中射血無力組EF 的均值和標(biāo)準(zhǔn)差為（0.377±0.093）%，射血正常組為（0.596±0.08）%。然后根據(jù)兩組受試者收縮時(shí)期的左心室位移場(chǎng)導(dǎo)出左心室在徑向、圓周方向上的最大累積位移、最高速度、最大累積應(yīng)變和最高應(yīng)變率。上述8 種描述左心室運(yùn)動(dòng)和形變特征的參數(shù)均由全局位移場(chǎng)計(jì)算得到，不再按照美國(guó)心臟協(xié)會(huì)（American Heart Association，AHA）的左心肌分段標(biāo)準(zhǔn)計(jì)算［19］。表1和表2 分別為兩組受試者的全局峰值位移和速度，以及全局峰值應(yīng)變和應(yīng)變率。

Table 1 Comparison of kinematic parameters between normal EF group and weak EF group表1 射血正常組與無力組受試者運(yùn)動(dòng)參數(shù)比較 （）

Table 1 Comparison of kinematic parameters between normal EF group and weak EF group表1 射血正常組與無力組受試者運(yùn)動(dòng)參數(shù)比較 （）

Table 2 Comparison of deformation parameters between normal EF group and weak EF group表2 射血正常組與無力組受試者形變參數(shù)比較（）

Table 2 Comparison of deformation parameters between normal EF group and weak EF group表2 射血正常組與無力組受試者形變參數(shù)比較（）

從表1、表2 中可以看出，射血正常組受試者運(yùn)動(dòng)和形變參數(shù)的絕對(duì)值均大于射血無力組，其中差距最明顯的為圓周應(yīng)變率（射血正常組為無力組的2 倍），其次為圓周應(yīng)變（射血正常組為無力組的1.78 倍），射血正常組受試者其余各項(xiàng)參數(shù)均為無力組的1.5 倍左右。射血正常組EF 為無力組的1.58 倍左右，這說明射血分?jǐn)?shù)越低，心力衰竭程度越高，從而導(dǎo)致左心室在徑向和圓周方向的總體收縮幅度越小，總體收縮力度越弱。此外，左心室EF 與其運(yùn)動(dòng)、形變參數(shù)存在一定相關(guān)性，具體如表3 所示。

Table 3 The correlation between EF of left ventricle and parameters of motion and deformation表3 左心室EF 與運(yùn)動(dòng)、形變參數(shù)的相關(guān)程度

由表3 可知，EF 與左心室在徑向、圓周方向上的運(yùn)動(dòng)和形變參數(shù)均具有較強(qiáng)相關(guān)性，其中相關(guān)程度最高的為徑向位移（相關(guān)系數(shù)為0.878），說明EF 受左心室徑向位移影響較大。這是由于徑向位移描述的是左心肌在收縮階段朝左心室中心運(yùn)動(dòng)的幅度，徑向位移越大，左心肌朝室中心擠壓得越厲害，左心室在收縮末期的體積越小，EF 也就越高。EF 正常受試者的左心室也可能存在運(yùn)動(dòng)異常情況，但本文計(jì)算的是整體位移場(chǎng)的運(yùn)動(dòng)與形變參數(shù)，未考慮局部節(jié)段的左心室運(yùn)動(dòng)。本文還對(duì)射血正常組和無力組受試者的各項(xiàng)運(yùn)動(dòng)與形變參數(shù)進(jìn)行了差異性分析，以驗(yàn)證建立的方法能否判斷左心室運(yùn)動(dòng)正常與否，具體如表4 所示。

Table 4 Difference of motion and deformation parameters between normal EF group and weak EF group表4 射血正常組與無力組受試者運(yùn)動(dòng)、形變參數(shù)差異程度

由表4 可知，兩組受試者各項(xiàng)參數(shù)均具有非常顯著的差異，說明本文方法求得的左心室位移場(chǎng)參數(shù)較為合理準(zhǔn)確，能在一定程度上判斷左心室是否正常運(yùn)動(dòng)。此外，為驗(yàn)證該方法的臨床表現(xiàn)，將實(shí)驗(yàn)結(jié)果與CVI 軟件進(jìn)行同種參數(shù)相關(guān)性比較，該軟件較為廣泛地應(yīng)用于臨床心臟功能評(píng)估［20］。表5 為本文方法與CVI 軟件各項(xiàng)參數(shù)的相關(guān)系數(shù)。

Table 5 Correlation analysis between the method and CVI software表5 本文方法與CVI 軟件的相關(guān)性分析

由表5 可知，除圓周方向上的位移和速度外，本文實(shí)驗(yàn)結(jié)果與CVI 軟件數(shù)據(jù)存在較高相關(guān)性，其中相關(guān)系數(shù)最高的為徑向位移（0.814）。不同的運(yùn)動(dòng)追蹤方法會(huì)因追蹤策略或參數(shù)求解方式不同而存在差異，例如CVI 軟件導(dǎo)出的圓周位移單位為弧度，而本文方法求取的圓周位移單位為毫米。進(jìn)一步比較CVI 軟件導(dǎo)出的射血正常組和無力組受試者圓周位移和速度的差異性，發(fā)現(xiàn)其差異系數(shù)分別為0.167 和0.029，不存在非常顯著的差異。

3 結(jié)語

量化左心室射血運(yùn)動(dòng)能為臨床診斷心血管疾病提供依據(jù)。本文首先利用三角網(wǎng)格生成左心室表面模型，形成輪廓點(diǎn)的連通網(wǎng)絡(luò)關(guān)系；然后擬合輪廓點(diǎn)局部曲面，計(jì)算輪廓點(diǎn)的局部形狀屬性；最后通過映射準(zhǔn)則恢復(fù)左心室位移場(chǎng)，并對(duì)初始位移場(chǎng)進(jìn)行平滑優(yōu)化。根據(jù)心臟電影磁共振短軸圖像數(shù)據(jù)，將101 例受試者分為射血正常組和射血無力組，比較了兩組受試者在徑向、圓周方向上的位移、速度、應(yīng)變和應(yīng)變率，結(jié)果表明本文建立的方法能有效區(qū)分左心室是否運(yùn)動(dòng)正常。對(duì)本文實(shí)驗(yàn)結(jié)果與CVI 心臟分析軟件數(shù)據(jù)進(jìn)行相關(guān)性分析，發(fā)現(xiàn)二者存在較高相關(guān)性，或可為臨床診療心血管疾病提供幫助。然而，由于難以找到合適且嚴(yán)格的標(biāo)準(zhǔn)描述位移場(chǎng)細(xì)節(jié)，對(duì)左心室位移場(chǎng)的驗(yàn)證仍存在一定難度和未知性［21］。