林劍圣，王麗嘉

結(jié)合生物力學(xué)模型重建左心室位移場

林劍圣，王麗嘉*

上海理工大學(xué) 醫(yī)療器械與食品學(xué)院，上海 200093

左心室心肌最為發(fā)達(dá)，心肌收縮產(chǎn)生的高壓將動脈血泵入全身，集中體現(xiàn)了心臟的泵血能力．定量分析左心室收縮運(yùn)動是診斷心血管疾?。ㄈ缧募」Ｋ溃┑闹匾緩剑疚牟捎妹枋鲎笮氖倚募〔馁|(zhì)的生物力學(xué)模型重建左心室位移場．該力學(xué)模型作為插值項(xiàng)，與心臟電影磁共振圖像的觀測位移場共同納入貝葉斯估計框架，并采用有限元法求解位移場方程．實(shí)驗(yàn)比較了左心室射血無力組（46例）與正常組（55例）的左心室功能參數(shù)，發(fā)現(xiàn)兩組在徑向和圓周方向的位移、速度、應(yīng)變和應(yīng)變率都具有非常顯著的差異（<0.001），這證明本文方法能夠有效區(qū)別左心室運(yùn)動正常與否．實(shí)驗(yàn)結(jié)果還與CVI軟件測量的左心室功能參數(shù)具有較高的相關(guān)性，說明本文方法有望輔助心血管疾病的臨床診斷．

左心室；位移場；生物力學(xué)模型；心臟電影磁共振圖像

引言

在世界范圍內(nèi)，尤其是低收入國家，心血管疾?。╟ardiovascular disease，CVD）仍然是導(dǎo)致死亡的主要原因[1]．在過去20年里，CVD也一直是我國城市和農(nóng)村地區(qū)人口死亡的主要原因之一[2]．如何有效預(yù)防和治療CVD刻不容緩．左心室（left ventricle，LV）是心臟泵血功能最重要的腔室，其心肌壁最厚、心肌收縮力最強(qiáng)．定量評估LV收縮功能，能夠輔助診斷心血管疾病，如缺血性心臟?。?/p>

臨床檢查心血管病的影像學(xué)手段包括心血管造影（angiocardiography）、超聲心動圖（echocardiography）、計算機(jī)斷層掃描（computed tomography，CT）、磁共振成像（magnetic resonance imaging，MRI）等．其中，心臟MRI除具有任意成像平面、心肌血流對比度高等獨(dú)特優(yōu)勢外，還能提供精確的心臟結(jié)構(gòu)和功能信息[3]．例如，心臟電影磁共振成像（cardiac cine MRI，CCMRI）因其高空間分辨率被廣泛用于射血分?jǐn)?shù)（ejection fraction，EF）、腔室體積、心肌質(zhì)量等的臨床測量．它還是無創(chuàng)量化心臟整體功能的金標(biāo)準(zhǔn)[4]．CCMRI的時間分辨率也足夠高，可以觀察心臟周期運(yùn)動．CCMRI可提供心臟解剖結(jié)構(gòu)，而標(biāo)記磁共振成像（tagged MRI，TMRI）[5]則是心臟MRI的功能成像類型．TMRI圖像的網(wǎng)格變化直觀記錄心肌形變，通常被視為二維層面評估心肌存活質(zhì)量的金標(biāo)準(zhǔn)．Chen等[6]利用Gabor濾波器檢測標(biāo)記網(wǎng)格點(diǎn)，運(yùn)用魯棒點(diǎn)匹配策略追蹤網(wǎng)格點(diǎn)運(yùn)動．Li等[7]利用圖割算法對Gabor濾波器檢測到的網(wǎng)格點(diǎn)分類，運(yùn)用有限差分法直接從分類網(wǎng)格點(diǎn)中計算雙心室三維形變．但TMRI圖像分辨率較低，而且網(wǎng)格線會隨著時間增長逐漸暗淡，影響運(yùn)動末期的參數(shù)提取．Yu等[8]提出結(jié)合稀疏先驗(yàn)約束的形變模型，該模型在TMRI圖像質(zhì)量較差、噪聲較多的情況下，能夠較好處理異常值和嚴(yán)重錯誤．其他用于心室功能分析的MRI技術(shù)包括相位對比、位移編碼和應(yīng)變編碼MRI，它們共同的缺點(diǎn)是在檢查時須獲得專門的圖像數(shù)據(jù)集，增加了檢查時間，且不能對現(xiàn)有的心臟磁共振數(shù)據(jù)集進(jìn)行回顧性應(yīng)用[9]．

從心臟磁共振圖像中評估心臟功能的前提一般是追蹤心臟運(yùn)動以及恢復(fù)位移場．傳統(tǒng)的運(yùn)動追蹤方法包括基于輪廓特征的方法、形變配準(zhǔn)和光流法等．Shi等[10]提倡把有實(shí)際物理意義的連續(xù)介質(zhì)力學(xué)原理納入到心臟運(yùn)動估計框架中，一方面基于輪廓點(diǎn)局部曲面在較小時間間隔發(fā)生微小形變的假設(shè)估計輪廓點(diǎn)最佳映射；另一方面把相位對比MRI融入CCMRI圖像中，作為LV室壁中運(yùn)動的補(bǔ)充信息，實(shí)現(xiàn)更加健壯完善的追蹤結(jié)果．Remme等[11]建立平面有限元網(wǎng)格模擬LV室壁的幾何變化：首先追蹤標(biāo)記的輪廓點(diǎn)運(yùn)動；然后將運(yùn)動結(jié)果與網(wǎng)格擬合，恢復(fù)LV三維形變；最后從TMRI圖像中導(dǎo)出運(yùn)動分布模型對原始形變進(jìn)行濾波，得到更加精確的結(jié)果．Veress等[12]提出超彈性翹曲（hyperelastic warping，HW）方法，利用增廣拉格朗日技術(shù)實(shí)現(xiàn)圖像能量硬約束最小化，同時利用超彈性應(yīng)變能量對圖像配準(zhǔn)實(shí)施軟約束．后者用來描述LV真實(shí)的本構(gòu)模型，該模型結(jié)合有限元離散化確保形變是微分同胚的，有利于追蹤心肌大范圍形變與旋轉(zhuǎn)運(yùn)動．心肌的近似不可壓縮性表現(xiàn)為心肌在有血流灌注情況下的體積變化不超過原體積的4%．Bistoquet等[13]將這一特性作為相鄰幀結(jié)點(diǎn)匹配應(yīng)當(dāng)遵循的約束條件，同時結(jié)合歸一化互信息配準(zhǔn)法在LV中間面曲線坐標(biāo)系中恢復(fù)LV三維周期運(yùn)動．Liu等[14]利用狀態(tài)空間分析技術(shù)構(gòu)建和整合心肌系統(tǒng)動力學(xué)、圖像觀測數(shù)據(jù)、以及處理和測量噪聲干擾．有實(shí)際物理意義的力學(xué)約束不僅作為心肌形變先驗(yàn)，而且插值和濾波圖像觀測數(shù)據(jù)，另外利用攜帶時空約束的統(tǒng)計濾波理論促進(jìn)合并多幀信息，以產(chǎn)生最優(yōu)的心臟運(yùn)動估計．心臟功能分析大多源于單一的功能或者結(jié)構(gòu)圖像，為了更好地表征特異的心臟生理病理，Wong等[15]利用心臟生理組模型有效融合差異較大的不同圖像源信息，這些圖像源包括CCMRI圖像的觀測位

移場和體表電位圖提供的心臟電生理時空特征，兩者互補(bǔ)可得到更加可靠準(zhǔn)確的心臟功能分析結(jié)果．

目前，基于深度學(xué)習(xí)的運(yùn)動追蹤存在挑戰(zhàn)，因?yàn)樾枰愵悢?shù)據(jù)訓(xùn)練模型泛化能力，而且真實(shí)的LV位移場幾乎不可能獲?。甉iao等[16]等比較了傳統(tǒng)配準(zhǔn)和卷積神經(jīng)網(wǎng)絡(luò)追蹤心臟運(yùn)動的性能，發(fā)現(xiàn)配準(zhǔn)方法性能更為穩(wěn)定，但耗時更長；而神經(jīng)網(wǎng)絡(luò)依賴訓(xùn)練數(shù)據(jù)，例如從健康志愿者隊(duì)列訓(xùn)練得到的模型用于追蹤病患心臟運(yùn)動，相比傳統(tǒng)配準(zhǔn)，其性能會顯著降低．Tang等[17]構(gòu)建密集連接Siamese卷積網(wǎng)絡(luò)，該網(wǎng)絡(luò)基于L2代價函數(shù)學(xué)習(xí)匹配最相似的LV局部面片特征．Zhang等[18]利用循環(huán)神經(jīng)網(wǎng)絡(luò)提取LV局部運(yùn)動特征，并采用全連接網(wǎng)絡(luò)識別心肌梗死區(qū)域．循環(huán)神經(jīng)網(wǎng)絡(luò)的訓(xùn)練和驗(yàn)證金標(biāo)準(zhǔn)是延遲釓增強(qiáng)MRI確定的心肌梗死區(qū)域．Zheng等[19]提出基于半監(jiān)督學(xué)習(xí)的心臟病理分類模型，該模型可根據(jù)心肌節(jié)段半徑和厚度特征區(qū)分不同類型的心血管疾?。鳰antilla等[20]采用字典學(xué)習(xí)區(qū)分LV運(yùn)動正常與否，判斷的關(guān)鍵特征是局部徑向時空輪廓．

從圖像中觀測的位移場因噪聲干擾顯得稀疏雜亂．為了恢復(fù)密集規(guī)律的位移場，有必要加入正則模型規(guī)范位移場秩序，例如提供平滑約束的自由形式形變（free form deformation，F(xiàn)FD）[21]和統(tǒng)計LV形變特征的主動形狀模型（active shape model，ASM）[22]．但本文作者認(rèn)為單憑圖像特征獲取的心臟運(yùn)動參數(shù)是不太可靠的，這是由CCMRI攜帶運(yùn)動信息稀疏所致．深度學(xué)習(xí)在醫(yī)學(xué)圖像分割領(lǐng)域展示了強(qiáng)大的性能[23,24]，但應(yīng)用于心臟運(yùn)動追蹤時存在局限性，因?yàn)闆]有一個統(tǒng)一嚴(yán)格的學(xué)習(xí)標(biāo)簽供機(jī)器訓(xùn)練[16-18]．多圖像源數(shù)據(jù)的互補(bǔ)[10,11,15]有利于實(shí)現(xiàn)更加準(zhǔn)確完善的運(yùn)動追蹤結(jié)果，但臨床更加傾向于常規(guī)的CCMRI檢查，CCMRI圖像精細(xì)的解剖結(jié)構(gòu)更易被患者和醫(yī)生采納，TMRI[6-8]更多地被用作科研序列．醫(yī)生可能更在意的是描述心肌局部運(yùn)動和形變的各種功能參數(shù)，以便輔助診斷心血管疾病，而非直接的分類結(jié)果[19,20]，這些參數(shù)無法由人眼視覺和臨床經(jīng)驗(yàn)直接得到．

本文采用描述LV心肌材質(zhì)屬性的生物力學(xué)模型（biomechanical model，BM）[25]，這與文獻(xiàn)[12-14]的研究思路類似，希望重建的LV位移場更符合實(shí)際的LV運(yùn)動特征．具體做法是將生物力學(xué)模型與從心臟電影磁共振圖像中觀測的位移場共同納入貝葉斯估計框架，然后采用有限元法求解力學(xué)模型引導(dǎo)下的優(yōu)化位移場．本文首先詳細(xì)介紹LV位移場重建流程，然后利用實(shí)驗(yàn)驗(yàn)證了本文所提方法的有效性．

1 理論方法

本文結(jié)合描述LV心肌材質(zhì)屬性的生物力學(xué)模型和心臟電影磁共振圖像的觀測位移場信息，重建LV完整位移場，整個方法流程如圖1所示．首先根據(jù)觀測位移場中的輪廓點(diǎn)位移向量和置信度構(gòu)造觀測位移場的高斯噪聲模型．在應(yīng)變能量函數(shù)中引入左心肌材質(zhì)矩陣，以便將LV應(yīng)變和形變勢能聯(lián)系起來．然后將能量項(xiàng)和噪聲項(xiàng)納入貝葉斯估計（Bayesian estimation）框架．為了求解優(yōu)化位移場，先對LV六面體剖分和四面體細(xì)分；再根據(jù)有限元理論，以矩陣形式重寫位移場方程，以便求出方程的一階導(dǎo)函數(shù)及LV總體剛度矩陣．最后以一定的迭代次數(shù)優(yōu)化觀測位移場，根據(jù)最終得到的LV優(yōu)化位移場，計算LV功能參數(shù)．

1.1 高斯噪聲模型

LV觀測位移場并非實(shí)際位移場，它會因成像缺陷、算法誤差等干擾因素而存在高斯噪聲[26]．假設(shè)從圖像中獲取觀測位移場的概率呈零均值正態(tài)分布，可建立如(1)式所示的高斯噪聲模型（Gaussian noise model）．

圖1 左心室位移場重建流程

由于本文的追蹤方法基于輪廓特征，因此精確的LV分割結(jié)果非常重要．心血管成像（circle cardiovascular imaging，CVI）軟件提供了穩(wěn)定且實(shí)時的全自動LV分割．但本文定義的LV輪廓是平滑非凹的，但CVI分割的原始輪廓是較為粗糙的，須進(jìn)一步對輪廓進(jìn)行閉合檢測、細(xì)化、平滑等處理．盡管CVI分割LV為本文科研工作提供巨大便利，但處理大量圖像難免會產(chǎn)生過分割、欠分割和漏分割個例．本文手動修復(fù)了錯誤分割個例，避免分割錯誤對功能參數(shù)估計的影響．除LV分割流程外，我們的三維建模、運(yùn)動追蹤以及應(yīng)變分析框架是完全自動穩(wěn)定的．有些不需要分割直接估計LV運(yùn)動參數(shù)的方法可能忽略了CCMRI圖像的一個特性，那就是室中層面的圖像幀存在較多的乳頭肌和小梁肌，干擾組織的突兀變化可能會導(dǎo)致跟蹤結(jié)果在解剖學(xué)上的不現(xiàn)實(shí)，這也是本文要求分割輪廓盡可能是凸的原因．此外，如果需要評估不同心肌節(jié)段的功能參數(shù)，首要條件是分割LV，然后依據(jù)美國心臟協(xié)會制定的左心肌分段標(biāo)準(zhǔn)再次細(xì)分LV．圖2(b)表面模型生成包括輪廓層間插值，以彌補(bǔ)層間分辨率小于層內(nèi)分辨率的差距，以及用三角面片表達(dá)輪廓點(diǎn)局部特征，該特征作為運(yùn)動追蹤的標(biāo)記，如圖2(c)所示．根據(jù)文獻(xiàn)[10]提供的薄板彈性能量在瞬時變化微小的假設(shè)完成前后幀輪廓點(diǎn)的匹配，最終得到LV觀測位移場，如圖2(d)所示．

1.2 應(yīng)變能量函數(shù)

物體形變產(chǎn)生的內(nèi)能變化可采用基于連續(xù)介質(zhì)力學(xué)的應(yīng)變能量函數(shù)表示[28]．該函數(shù)不考慮剛體運(yùn)動如平移、旋轉(zhuǎn)等對物體形變的影響，如(7)式所示．

其中，D是物體材質(zhì)屬性矩陣，它將物體的應(yīng)變E同物體內(nèi)能聯(lián)系起來．應(yīng)變E是由三個正應(yīng)變和三個剪應(yīng)變組成的列向量，反映了物體的形變程度和類型[29]．本文采用的生物力學(xué)模型是將LV心肌看成橫向各向同性材質(zhì)，考慮了肌纖維方向的優(yōu)先剛度，其材質(zhì)矩陣D的逆矩陣如(8)式所示．

接著，需要把應(yīng)變能量函數(shù)轉(zhuǎn)化為貝葉斯估計中的先驗(yàn)概率．物體在某一點(diǎn)處是否發(fā)生形變是通過與它相鄰質(zhì)點(diǎn)的相對位置變化體現(xiàn)的．而應(yīng)變能量函數(shù)通過應(yīng)變E反映了某點(diǎn)處的形變勢能，即某點(diǎn)應(yīng)變能量由該點(diǎn)與鄰域點(diǎn)在物體運(yùn)動前后的位置關(guān)系決定．概率密度函數(shù)中類似的是馬爾科夫隨機(jī)場（Markov random field），場中某點(diǎn)的概率值由鄰域點(diǎn)決定[30]．因此，應(yīng)變能量函數(shù)的概率形式可寫為：

其中，為歸一化常量．

1.3 貝葉斯估計

以貝葉斯估計形式將圖像觀測位移場與生物力學(xué)模型結(jié)合起來，重建LV完整位移場．即把觀測位移場的高斯噪聲模型和引入LV心肌材質(zhì)屬性的應(yīng)變能量函數(shù)共同納入貝葉斯估計框架中，如(10)式所示．

接著，將(1)式和(9)式代入(11)式，省去符號，忽略常數(shù)項(xiàng)和不包含的項(xiàng)，可得：

最后，采用有限元法求解位移場方程之前，需要先對LV進(jìn)行有限元剖分．

1.4 LV有限元剖分

有限元剖分可以逼近LV復(fù)雜的幾何形狀．本文先采用對稱最近鄰法連接同一層面的LV輪廓點(diǎn)（包括心內(nèi)膜、心肌和心外膜輪廓），對稱最近鄰點(diǎn)對的關(guān)系表達(dá)式為：

其中，是歐式距離算子，和表示不同輪廓類型的點(diǎn)索引．如果點(diǎn)的最近鄰點(diǎn)為，同時點(diǎn)的最近鄰點(diǎn)為，那么和的關(guān)系是對稱最近鄰．接著，以自下而上的方式（從心尖到心底的順序）進(jìn)行同一輪廓類型的點(diǎn)采樣，采樣方案如下：

其中，和分別表示剩余輪廓點(diǎn)數(shù)和剩余采樣輪廓點(diǎn)數(shù)，int為取整運(yùn)算符，和分別表示當(dāng)前輪廓點(diǎn)索引和采樣輪廓點(diǎn)索引．最后，配對不同輪廓類型的采樣輪廓點(diǎn)，配對方案如下：

其中，ds是以p和為起始點(diǎn)對的歐式距離和算子，Emp表示使得所有配對輪廓點(diǎn)歐式距離和最小的起始點(diǎn)對．上述步驟分別完成了不同輪廓類型的點(diǎn)連接，同一輪廓類型的點(diǎn)采樣以及不同輪廓類型的點(diǎn)配對，最終得到緊密連接的LV六面體網(wǎng)格，如圖3所示．

本文將六面體網(wǎng)格的8個結(jié)點(diǎn)進(jìn)行編號，如圖4(a)所示．四面體網(wǎng)格是在六面體網(wǎng)格的基礎(chǔ)上繼續(xù)細(xì)分，以六面體網(wǎng)格結(jié)點(diǎn)編號為參考，四面體網(wǎng)格結(jié)點(diǎn)編號方式有6種，分別為（1-2-5-8）、（2-5-7-8）、（2-5-6-7）、（1-3-4-8）、（1-3-7-8）、（1-2-3-7）．根據(jù)上述編號方式，可以將一個六面體網(wǎng)格繼續(xù)剖分為6個四面體網(wǎng)格，如圖4(b)和4(c)所示．LV剖分過程中可能會出現(xiàn)病態(tài)四面體、退化四面體，這是不良結(jié)構(gòu)，影響形變參數(shù)的提取，而且剖分過程沒有考慮到患者LV的特異性．文獻(xiàn)[31]提供了一種無網(wǎng)格法可用于解決心臟大形變，材質(zhì)不均勻等困難，這將是本文應(yīng)變分析模塊的改進(jìn)方向．

圖4 六面體網(wǎng)格編號與四面體網(wǎng)格剖分

1.5 有限元法解方程

有限元法中的離散化步驟是利用形狀函數(shù)將結(jié)點(diǎn)位移和單元的內(nèi)位移聯(lián)系起來．四面體單元內(nèi)一點(diǎn)的位移可由4個結(jié)點(diǎn)的形狀函數(shù)和結(jié)點(diǎn)位移來表示．

根據(jù)有限元理論，可從(18)式中獲取總體剛度矩陣K：

其中，C是總體置信度矩陣，為一對角矩陣．將(20)式和(21)式代入(12)式中，得到矩陣形式的貝葉斯估計方程，如(22)式所示．

其中任一子矩陣為：

其中，r和s是四面體單元結(jié)點(diǎn)索引之一（i、j、k、m）．和是對應(yīng)結(jié)點(diǎn)的應(yīng)變矩陣，矩陣中的每項(xiàng)元素都是由結(jié)點(diǎn)坐標(biāo)決定的常數(shù)[32]．將每個按照結(jié)點(diǎn)索引順序投放，最終形成總體剛度矩陣K．把總體剛度矩陣K、總體置信度矩陣C和觀測位移場代入(23)式，并以適當(dāng)?shù)拇螖?shù)迭代求解，將從觀測位移場中得到LV優(yōu)化位移場，如圖5所示．圖5(a)和5(c)分別是收縮和舒張時期的觀測位移場，圖5(b)和5(d)分別是收縮和舒張時期對應(yīng)的優(yōu)化位移場．仔細(xì)觀察圖5可以發(fā)現(xiàn)，LV優(yōu)化位移場比原始觀測位移場更加密集有序．

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

本文使用的101例心臟電影磁共振圖像數(shù)據(jù)由46例LV射血無力患者和55例LV射血正?；颊呓M成，射血無力表現(xiàn)為左心室射血分?jǐn)?shù)（LVEF）低于50%．射血正常和射血無力的都是心血管疾病患者（主要是心肌肥厚、心力衰竭、心室重構(gòu)患者）．患者性別包括35例女性和66例男性，年齡處于14～88歲，心率范圍為46～125次/min．磁共振圖像在GE 1.5T掃描儀上通過穩(wěn)態(tài)自由進(jìn)動（steady state free precession，SSFP）成像序列獲取，具體成像參數(shù)如下：圖像矩陣大小為256×256，視野為360 mm ×360 mm，重復(fù)時間為3.5 ms，回波時間為1.4 ms，層厚為6～7 mm，層間距為7～10 mm．?dāng)?shù)據(jù)來源于康奈爾醫(yī)學(xué)院，患者已簽訂知情同意書，通過了上海交通大學(xué)醫(yī)學(xué)院附屬上海兒童醫(yī)學(xué)中心的倫理審批要求．后續(xù)數(shù)據(jù)處理在MATLAB2020A上完成，LVEF首先是根據(jù)生成的LV三維模型，通過convhull函數(shù)估計舒張末期與收縮末期體積，然后依據(jù)EF公式計算得到．

3 結(jié)果與討論

3.1 LV功能參數(shù)

本文從重建的LV位移場中導(dǎo)出了LV功能參數(shù)，包括徑向位移（radial displacement，RD）和速率（radial velocity，RV）、圓周位移（circumferential displacement，CD）和速率（circumferential velocity，CV），徑向應(yīng)變（radial strain，RS）和應(yīng)變率（radial strain rate，RSR）、圓周應(yīng)變（circumferential strain，CS）和應(yīng)變率（circumferential strain rate，CSR）．RD和RS描述LV的收縮和舒張程度，CD和CS刻畫LV的扭轉(zhuǎn)和解旋程度．LV在一個心動周期的全局累積RD和CD曲線如圖6左上所示，全局累積RS和CS曲線如圖6右上所示．全局累積曲線表示LV從舒張末期開始，逐步收縮和扭轉(zhuǎn)，一直到收縮末期的收縮和扭轉(zhuǎn)峰值狀態(tài)，再逐步舒張和解旋，最終重新回到舒張末期的原始狀態(tài)．圖6左下的RV和CV曲線表示每幀時間間隔內(nèi)，LV在徑向和圓周方向上的整體運(yùn)動快慢，圖6右下的RSR和CSR曲線表示每幀時間間隔內(nèi)，LV在徑向和圓周方向上的整體形變快慢．

圖6 左心室功能參數(shù)曲線

為了驗(yàn)證本文方法能夠有效區(qū)分LV收縮運(yùn)動正常與否，本節(jié)比較了LV射血正常組和無力組的LV各功能參數(shù)的差異，包括徑向和圓周方向上的峰值位移和應(yīng)變，以及收縮時期的最大速率和應(yīng)變率（參照圖6各分圖收縮期全局峰值點(diǎn)），如表1所示．

表1 射血正常組與無力組的左心室功能參數(shù)比較

由表1可知，射血正常組的各項(xiàng)LV功能參數(shù)的絕對值要比射血無力組更大，在一定程度上說明LV的EF和其整體收縮能力存在正相關(guān)：EF越高，其整體收縮能力也就越強(qiáng)．但也要注意到，由于本文未局部測量LV功能參數(shù)，因此射血正?；颊叩腖V功能參數(shù)也可能出現(xiàn)偏低的情況．兩組各項(xiàng)參數(shù)的值均小于0.001，說明兩組的LV功能參數(shù)存在非常顯著的差異，證明了本文方法能夠有效區(qū)分LV整體收縮運(yùn)動正常與否．

3.2 本文方法與CVI軟件比較

為了進(jìn)一步驗(yàn)證本文所提方法的臨床實(shí)用性，將本文方法測量的LV功能參數(shù)與心血管成像軟件CVI進(jìn)行相關(guān)性分析，該軟件較為廣泛地應(yīng)用于臨床心臟功能分析以及心血管疾病的科學(xué)研究．圖7顯示了本文方法測量的101例收縮時期全局峰值參數(shù)包括徑向應(yīng)變、應(yīng)變率、位移、速度以及圓周應(yīng)變、應(yīng)變率與CVI軟件測量結(jié)果的相關(guān)性，對應(yīng)的線性回歸方程均在圖中標(biāo)出．RS、RSR、RD、RV、CS、CSR的相關(guān)系數(shù)分別為0.718、0.726、0.668、0.765、0.708和0.792，其中CSR的相關(guān)性最高．高相關(guān)性說明了本文方法與CVI軟件的匹配程度較高（0.6～0.8時為強(qiáng)相關(guān)）．由于兩者測量的圓周位移和速率單位不一致，本文方法為毫米，CVI為弧度，所以這兩項(xiàng)參數(shù)沒有相關(guān)性．此外，CVI測量的LV射血正常組和無力組的圓周位移和速率不存在非常顯著的差異（詳見附件表S1，掃描文章首頁OSID二維碼，或在網(wǎng)頁版獲?。?/p>

圖7 使用本文方法與CVI軟件獲取的左心室功能參數(shù)的相關(guān)性

3.3 不同LVEF患者的LV整體與局部功能比較

分別隨機(jī)選取2例臨床表現(xiàn)為LV射血無力和射血正常對象，比較他們的RD、RS、CD、CS在整個心動周期的演化趨勢，如圖8所示．累積應(yīng)變和相同方向的累積位移是關(guān)系密切的，LVEF也和累積峰值存在一定的相關(guān)性，如果不考慮射血分?jǐn)?shù)保留型心衰疾病的話．LVEF越高，功能參數(shù)的峰值越大，電影幀表現(xiàn)為心肌增厚量更多．

圖8 具有不同LVEF研究對象的左心室功能參數(shù)曲線

由上述分析可知，LVEF是LV整體功能指標(biāo)，它和全局累積徑向位移、應(yīng)變存在正相關(guān)，和全局累積圓周位移、應(yīng)變存在負(fù)相關(guān)．但是心肌可能存在局部運(yùn)動異常，這可能無法從全局功能參數(shù)和LVEF得到．繼續(xù)觀察1例左心室射血正常（LVEF>50%）與1例射血無力（LVEF<50%）患者在收縮時期的應(yīng)變云圖，如圖9所示．主應(yīng)變（顏色棒范圍值）是結(jié)點(diǎn)徑向和圓周應(yīng)變平方和的正平方根，用來衡量結(jié)點(diǎn)局部區(qū)域總形變程度，平均主應(yīng)變是所有結(jié)點(diǎn)主應(yīng)變的和除以所有結(jié)點(diǎn)數(shù)量．由圖9可知，從局部層面看，收縮階段的不同幀LV射血正常組大部分心肌節(jié)段的形變量明顯大于射血無力患者，而且存在較多大形變（紅色部分），該形變量是指當(dāng)前幀到下幀的瞬時變化．此外從整體上看，LV射血正常患者組的結(jié)點(diǎn)平均主應(yīng)變均大于同時刻的射血無力患者．射血正常組的LV形態(tài)變化較大，心肌增厚量多，而射血無力患者的LV形態(tài)變化小，心肌厚度幾乎不變．

圖9 具有不同LVEF研究對象的左心室應(yīng)變云圖

3.4 不同約束條件下的LV功能參數(shù)比較

LV位移場在不同約束條件的位移大小和方向會發(fā)生變化，其功能參數(shù)也會相應(yīng)變化．101例患者原始位移場、局域平滑位移場和生物力學(xué)模型指導(dǎo)位移場的LV功能參數(shù)差異如表2所示．分別對比無約束和模型約束的RS和CS絕對值，可以發(fā)現(xiàn)無約束的RS比模型約束的大，CS比模型約束的小，而兩者總應(yīng)變量差不多．因此生物力學(xué)模型提供了LV的旋轉(zhuǎn)動力，將更多的形變量分配到周向，增加更多周向位移．平滑約束會稍微降低總形變量．根據(jù)進(jìn)一步實(shí)驗(yàn)得出不同約束條件下的射血無力和正常組的LV功能參數(shù)絕大部分存在非常顯著的差異（詳見附件表S1），說明不同類型位移場都可用于區(qū)分LV正常與否，證明本文運(yùn)動追蹤方法的可靠性．而生物力學(xué)模型更多地是增加LV旋轉(zhuǎn)運(yùn)動，使位移場演化更符合真實(shí)的LV運(yùn)動．因?yàn)檎鎸?shí)LV運(yùn)動是存在解旋和扭轉(zhuǎn)的，這點(diǎn)應(yīng)從圓周方向的參數(shù)看出．追蹤算法導(dǎo)出的初始位移場是較為雜亂無章的，力學(xué)模型在原始位移場的基礎(chǔ)上，增加了一些周向的分量．

表2 不同約束條件的左心室功能參數(shù)比較

4 結(jié)論

為了準(zhǔn)確提取心臟電影磁共振圖像中的LV運(yùn)動信息，輔助診斷心血管疾病，本文將描述LV心肌材質(zhì)屬性的生物力學(xué)模型作為心臟電影磁共振圖像觀測位移場的插值項(xiàng)，重建LV優(yōu)化位移場．具體做法是將兩者共同納入貝葉斯估計框架，然后采用有限元法解位移場方程，最后依據(jù)優(yōu)化位移場，導(dǎo)出LV功能參數(shù)．為了驗(yàn)證本文方法能夠有效區(qū)分LV運(yùn)動正常與否，本文比較了LV射血正常組和無力組的LV功能參數(shù)，發(fā)現(xiàn)兩組在徑向和圓周方向上的位移、速度、應(yīng)變和應(yīng)變率均具有非常顯著的差異．本文方法還與CVI軟件的測量結(jié)果具有較高的相關(guān)性，說明本文方法有望輔助臨床診療心血管疾?。?/p>

無

附件材料

表S1 不同約束條件下，本文方法和CVI軟件測量的射血正常和射血無力患者各LV功能參數(shù)的差異

[1] DEMARINIS S.Cancer overtakes cardiovascular disease as leading cause of death in wealthy nations[J].Explore, 2020, 16(1): 6-7.

[2] PANG Y, LYU J, YU C, et al.Risk factors for cardiovascular disease in the Chinese population: recent progress and implications[J].Glo Hea J, 2020, 4(3): 65-71.

[3] FRANGI A F, NIESSEN W J, VIERGEVER M A.Three-dimensional modeling for functional analysis of cardiac images: a review[J].IEEE Trans Med Imaging, 2001, 20(1): 2-5.

[4] WANG H, AMINI A A.Cardiac motion and deformation recovery from MRI: A Review[J].IEEE Trans Med Imaging, 2012, 31(2): 487-503.

[5] ZERHOUNI E A, PARISH D M, ROGERS W J, et al.Human heart: Tagging with MR imaging - A method for noninvasive assessment of myocardial motion[J].Radiology, 1988, 169(1): 59-63.

[6] CHEN T, WANG X, CHUNG S, et al.Automated 3D motion tracking using Gabor filter bank, robust point matching, and deformable models[J].IEEE Trans Med Imaging, 2009, 29(1): 1-11.

[7] LI M, GUPTA H, LLOYD S G, et al.A graph theoretic approach for computing 3D+ time biventricular cardiac strain from tagged MRI data[J].Med Image Anal, 2017, 35: 46-57.

[8] YU Y, ZHANG S T, LI K, et al.Deformable models with sparsity constraints for cardiac motion analysis[J].Med Image Anal, 2014, 18(6): 927-937.

[9] TSADOK Y, FRIEDMAN Z, HALUSKA B A, et al.Myocardial strain assessment by cine cardiac magnetic resonance imaging using non-rigid registration[J].Magn Reson Imaging, 2016, 34(4): 381-390.

[10] SHI P, SINUSAS A J, CONSTABLE R T, et al.Volumetric deformation analysis using mechanics-based data fusion: applications in cardiac motion recovery[J].Int J Comput Vision, 1999, 35(1):87-107.

[11] REMME E W, AUGENSTEIN K F, YOUNG A A, et al.Parameter distribution models for estimation of population based left ventricular deformation using sparse fiducial markers[J].IEEE Trans Med Imaging, 2005, 24(3): 381-388.

[12] VERESS A I, GULLBERG G T, WEISS J A.Measurement of strain in the left ventricle during diastole with cine-MRI and deformable image registration[J].J Biomech Eng, 2005, 127(7): 1195-1207.

[13] BISTOQUET A, SKRINJAR O M.Left ventricular deformation recovery from cine MRI using a 4d incompressible model[J].IEEE Trans Med Imaging, 2007, 26(9): 1136-1153.

[14] LIU H F, SHI P C.State-space analysis of cardiac motion with biomechanical constraints[J].IEEE T Image Process, 2007, 16(4): 901-917.

[15] WONG K C L, WANG L W, ZHANG H Y, et al.Physiological fusion of functional and structural images for cardiac deformation recovery[J].IEEE Trans Med Imaging, 2011, 30(4): 990-1000.

[16] QIAO M Y, WANG Y Y, GUO Y, et al.Temporally coherent cardiac motion tracking from cine MRI: traditional registration method and modern CNN method[J].Med Phys, 2020, 47(9): 4189-4198.

[17] TANG J Y, GAN Z Y, YANG X.Cardiac motion tracking in short-axis MRI using siamese convolution network[C]//2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM).IEEE, 2019: 865-870.

[18] ZHANG N, YANG G, GAO Z F, et al.Deep learning for diagnosis of chronic myocardial infarction on nonenhanced cardiac cine MRI[J].Radiology, 2019, 291(3): 606-617.

[19] ZHENG Q, DELINGETTE H, AYACHE N.Explainable cardiac pathology classification on cine MRI with motion characterization by semi-supervised learning of apparent flow[J].Med Image Anal, 2019, 56: 80-95.

[20] MANTILLA J J, PAREDES J L, BELLANGER J J, et al.Discriminative dictionary learning for local LV wall motion classification in cardiac MRI[J].Expert Syst Appl, 2019, 129: 286-295.

[21] SEDERBERG T W, PARRY S R.Free-form deformation of solid geometric models[J].Siggraph, 1986, 20(4): 151-160.

[22] COOTES T F, TAYLOR C J, COOPER D H, et al.Active shape models-their training and application[J].Comput Vis Image Underst, 1995, 61(1): 38-59.

[23] GONG J C, WANG Y, WANG Y J.A method for segmentation of glioma on multimodal magnetic resonance images based on wavelet fusion and deep learning[J].Chinese J Magn Reson, 2020, 37(2): 131-143.

宮進(jìn)昌, 王宇, 王遠(yuǎn)軍.結(jié)合小波融合和深度學(xué)習(xí)的腦膠質(zhì)瘤自動分割[J].波譜學(xué)雜志, 2020, 37(2): 131-143.

[24] LIU P, ZHONG Y M, WANG L J.Automatic segmentation of right ventricle in cine cardiac magnetic resonance image based on a dense and multi-scale U-net method[J].Chinese J Magn Reson, 2020, 37(4): 456-468.

劉鵬, 鐘玉敏, 王麗嘉.基于密集多尺度U-net網(wǎng)絡(luò)的電影心臟磁共振圖像右心室自動分割[J].波譜學(xué)雜志, 2020, 37(4): 456-468.

[25] GUCCIONE J M, MCCULLOCH A D.Finite element modeling of ventricular mechanics[M].Springer New York, 1991.

[26] PAPADEMETRIS X, SHI P C, DIONE D P, et al.Recovery of soft tissue object deformation from 3D image sequences using biomechanical models[C]// Proceedings of the 16th International Conference on Information Processing in Medical Imaging.Berlin, Heidelberg: Springer, 2000.

[27] SHI P C, SINUSAS A J.Point-tracked quantitative analysis of left ventricular surface motion from 3-D image sequences[J].IEEE Trans Med Imaging, 2000, 19(1): 36-50.

[28] PAPADEMETRIS, XENOPHON, SINUSAS, et al.Estimation of 3-D left ventricular deformation from medical images using biomechanical models[J].IEEE Trans Med Imaging, 2002, 21(7): 786-800.

[29] 楊桂通.彈性力學(xué)第2版[M].高等教育出版社, 2012.

[30] GEMAN S, GERMAN D.Stochastic relaxation, gibbs distributions and the Bayesian restoration of images[J].IEEE Trans Pattern Anal Mach Intell, 1984, 6: 721-741.

[31] ZHANG H Y, GAO Z, XU L, et al.A meshfree representation for cardiac medical image computing[J].IEEE J Transl Eng He, 2018, 6: 1-12.

[32] 李亞智.有限元法基礎(chǔ)與程序設(shè)計[M].科學(xué)出版社, 2004.

Reconstruction of Displacement Field of Left Ventricle Combined with Biomechanical Model

，WANG Li-jia

School of Medical Instrument and Food Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China

The left myocardium is the strongest cardiac muscle, representing the ability of heart to pump arterial blood throughout the body.Quantitative analysis of the contractive motion of left ventricle (LV) stands an important approach to diagnose cardiovascular diseases such as myocardial infarction.This paper utilized a biomechanical model describing the material of left myocardium to reconstruct the displacement field of LV.The mechanical model was used as an interpolation term, which was incorporated into the Bayesian estimation framework with the displacement field observed by cardiac cine magnetic resonance (CCMR) image, and the equation of displacement field was solved by finite element method.Functional parameters of LV were compared between the weak (46 cases) and normal (55 cases) ejection fraction groups of LV in the experiment.Significant differences in radial and circumferential displacement and velocity, strain and strain rate were observed between the two groups (<0.001), which proved that the method could effectively distinguish whether the motion of LV is normal or not.The experimental results are also highly correlated with the functional parameters of LV measured by CVI software, indicating that the method will facilitate the diagnosis of cardiovascular diseases.

left ventricle, displacement field, biomechanical model, cardiac cine magnetic resonance image

TP391

A

10.11938/cjmr20212898

2021-03-22;

2021-07-05

上海理工大學(xué)2021醫(yī)工交叉項(xiàng)目(10-21-308-412).

* Tel: 021-55271116; E-mail: lijiawangmri@163.com.