鄭南南，黃 鋼，林劍圣

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

0 引言

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

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

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

1 實驗方法

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

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

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

Fig.1 Framework of motion tracking algorithm of left ventricle圖1 左心室運動追蹤方法框架

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

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

1.2 左心室表面模型生成

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

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

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

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

1.3 特征追蹤

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

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

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

計算出系數(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）

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

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

將輪廓點局部曲面勢能映射到其鄰域范圍，得到輪廓的勢能分布情況，如圖5 所示。由圖5 可知，左心室輪廓的整體勢能較低，舒張期輪廓勢能比收縮期更低，說明舒張期輪廓的形變程度更小。

Fig.5 Potential energy distribution of contours diastole period（a）and systolic period（b）圖5 輪廓勢能分布舒張期（a）與收縮期（b）

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

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

1.4 位移場優(yōu)化

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

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

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

2 實驗結(jié)果與分析

2.1 實驗數(shù)據(jù)

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

2.2 實驗結(jié)果

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

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

Table 1 Comparison of kinematic parameters between normal EF group and weak EF group表1 射血正常組與無力組受試者運動參數(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 中可以看出，射血正常組受試者運動和形變參數(shù)的絕對值均大于射血無力組，其中差距最明顯的為圓周應變率（射血正常組為無力組的2 倍），其次為圓周應變（射血正常組為無力組的1.78 倍），射血正常組受試者其余各項參數(shù)均為無力組的1.5 倍左右。射血正常組EF 為無力組的1.58 倍左右，這說明射血分數(shù)越低，心力衰竭程度越高，從而導致左心室在徑向和圓周方向的總體收縮幅度越小，總體收縮力度越弱。此外，左心室EF 與其運動、形變參數(shù)存在一定相關性，具體如表3 所示。

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

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

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

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

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

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

3 結(jié)語

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