實時自由呼吸心臟成像是一種重復性好、采集時間短且不需要屏氣的磁共振成像方法。然而實時心臟成像進行左心室功能分析時, 通常需要依靠人眼找出舒張末期和收縮末期圖像, 這個過程非常費時費力。為了節省處理時間, 本文提出一種半自動定位舒張末期和收縮末期圖像的方法。該方法主要通過對呼吸運動信號進行擬合, 在擬合信號呼吸末期利用互相關系數確定舒張末期和收縮末期圖像。整個過程采集了10名健康志愿者實時心臟磁共振圖像, 并同標準屏氣穩態自由進動序列采集數據進行左心室功能分析比較。實驗結果表明, 所提方法能夠正確檢測出舒張末期和收縮末期圖像。未來, 這個方法可望在臨床上進行快速地左心室功能分析。
引用本文： 楊帆, 何艷, 張潔, 吳垠. 實時心臟磁共振成像左心室功能分析研究. 生物醫學工程學雜志, 2015, 32(6): 1279-1283. doi: 10.7507/1001-5515.20150227 復制
0 引言
心臟磁共振成像是一種精準的可重復用于左心室(left ventricle, LV)功能分析的方法[1]。目前，實時自由呼吸心臟成像可用來進行左心室的功能分析。與標準的屏氣穩態自由進動序列(steady-state free precession, SSFP)相比，實時成像能縮短數據采集時間且不需要屏氣，有助于不能屏氣的患者[2]。
然而，實時自由呼吸心臟成像左心室功能分析非常耗時，因為需要從多幅磁共振序列圖像中定位心臟舒張末期(end-diastole, ED)和收縮末期(end-systole, ES)圖像。并且由于呼吸運動，左心室的位置隨著橫膈膜的運動而發生變化，如果采用每層不同心臟位置的圖像進行數據分析，可能會導致左心室功能分析數據出錯, 其主要原因是橫膈膜運動導致每層心臟處于不同位置而產生形變。為了減少功能分析誤差，選擇相同心臟位置圖像，能有效解決這個問題。而呼氣末期相對吸氣階段，橫膈膜運動較平緩，因此可以使用左心室位于呼氣末期的ED和ES圖像分析心臟功能。目前，ED和ES圖像的選取主要采用人眼觀察的方法，從每層心臟磁共振序列中找出ED和ES圖像[3]。這種方法增加了后處理時間且非常費力，因此如果可自動定位ED和ES圖像，有利于節省數據處理時間，便于后續數據分析。當前在一個心動周期自動定位ED和ES圖像的方法[4]，需要受試者屏氣，而實時心臟磁共振成像中橫膈膜的運動導致每層心臟位置不同，因此無法直接采用該方法進行測定。
本文旨在提供一種能自動檢測實時心臟磁共振成像中ED和ES圖像的方法。由于心臟的位置隨著橫膈膜的運動而發生變化，左心室的運動軌跡可以通過呼吸運動信號進行擬合, 然后運用互相關系數檢測呼氣末期區域ED和ES圖像。利用本文方法僅需較少人工干預，即可精準檢測出ED和ES圖像，并用于左心室功能分析。
1 數據來源
本研究從中國科學院深圳先進技術研究院招募10名健康志愿者，年齡為(25±4)歲，無既往心臟病史，所有受試者均簽署知情同意書，并經過單位倫理審查，同意試驗。研究使用中國科學院深圳先進技術研究院3.0 T TIM TRIO西門子磁共振儀進行掃描。首先使用標準二維SSFP的具有心電圖(electrocardiogram，ECG)觸發和呼吸抑制的序列掃描心臟，10個短軸位層包含心臟從心尖到基底；實時自由呼吸的SSFP序列在短軸位對相同層位置進行掃描，并在時間方向運用K-L變換(Karhunen-Loeve Transform，KLT)濾波增加信噪比[5]。成像參數為：層厚8 mm，間隙2 mm，視野為340×287 mm2。對于標準SSFP心臟成像序列：重復時間/回波時間=3.4/1.5 ms，矩陣大小為256×216，并行采集(integrated parallel acquisition techniques, iPAT)加速因子為2，帶寬977 Hz/像素，時間分辨率40.7 ms。對于實時心臟成像序列：重復時間/回波時間=2.5/1.1 ms，矩陣大小為160×128，時間并行采集(temporal parallel acquisition techniques，TPAT)加速因子為4，帶寬1 488 Hz/像素，時間分辨率59.5 ms，并且掃描每層持續時間為5 s，包含84幅從呼氣末期到吸氣末期心臟圖像。
2 方法
2.1 左心室定位和分割
左心室定位首先需要確定心臟在磁共振圖像上的位置。以往定位心臟的方法主要針對屏氣心臟成像，并且假設心臟是唯一運動的器官[6-8]，而實時成像中，受呼吸運動影響，心臟不是唯一的運動器官，因此前述方法不能用于實時心臟成像圖像的心臟定位。文中采用人工定位的方法，對每層所有圖像求和，手動選取求和后圖像中左心室血池中心點或者離中心距離最近較亮的點做為種子點，如圖 1所示。圖 1左圖為1名志愿者心臟短軸位求和后圖像，右圖為種子點選取圖像，該點通過手動點擊獲取，并記錄坐標。以往左心室心內輪廓分割方法采用最優閾值，能有效實現左心室血池分割[9-10]。本研究同樣采用Otsu方法[11]計算每幅圖像分割閾值，將種子點作為每幅圖像區域增長的中心，相鄰點高于閾值加入種子點，低于閾值作為心內輪廓。通過該法能快速分割每層圖像左心室血池。左心室ED和ES血池分割結果如圖 2所示。
圖1
心臟磁共振短軸位求和圖像
Figure1.
Total image of short axis of cardiac magnetic resonance
圖2
左心室血池分割結果
Figure2.
Segmentation result of LV blood pool
2.2 橢圓擬合和呼吸運動擬合
為了降低乳突肌和肌小梁對于左心室中心檢測的影響，橢圓擬合之前先對分割血池進行凸包處理, 獲得最小多邊形邊界。如圖 3所示，左圖為左心室血池分割結果，右圖為運用凸包處理獲得的邊界。然后運用最小二乘法對左心室橢圓擬合，定位左心室中心。橢圓擬合圖如圖 4所示，其中藍色為擬合的橢圓，綠色為凸包邊界，紅色為左心室中心。
圖3
血池和凸包邊界
Figure3.
Blood pool and convex hull border
圖4
ED和ES橢圓擬合
Figure4.
Ellipse fitting of ED and ES
由于心臟位置隨著橫膈膜運動而變化，左心室在呼氣末期相對位置高于吸氣末期，可用心臟空間位置描述呼吸運動軌跡。然而由于心臟跳動，信號曲線不平滑，且考慮到心臟跳動頻率遠高于呼吸頻率，因此采用低通濾波器平滑曲線以獲得呼吸運動信號，漢寧窗降低信號波動。如圖 5所示，為1名志愿者呼吸運動擬合曲線，其中橫坐標表示每層心臟圖像序列，縱坐標表示左心室中心相對圖像邊界距離，心臟和呼吸運動信號用不同符號在圖中標記，AB區間為呼氣末期區間，包含一個心動周期。
圖5
呼吸運動擬合曲線
Figure5.
Fitting curve of respiratory motion
在心尖位置，由于左心室血池面積較小，初始種子點位置固定，呼吸運動可能導致收縮末期種子點位于血池之外，如圖 6所示。為了獲得呼吸運動信號，只對種子點位于血池內圖像進行橢圓擬合，記錄中心點位置，而血池外點排除，其中心位置通過周圍點插值獲得。最后通過多次信號擬合獲得呼吸運動信號。而在基底層，當心肌包圍左心室腔體面積超過50%時，進行橢圓擬合定位左心室中心，其它圖像的左心室中心位置同樣通過插值擬合獲得。基底層連續四個心臟周期的左心室ED和ES擬合圖像如圖 7所示，其中藍色為橢圓擬合曲線、紅點為左心室中心和綠色為凸包邊界。
圖6
左心室心尖收縮末期
Figure6.
ES of LV apex slice
圖7
基底層ED和ES擬合結果
Figure7.
Fittings results of ED and ES in basal slice
2.3 確定呼氣末期
定位ED和ES圖像之前需要確定呼氣末期位置。在本研究中，心跳范圍為60~100次/min，時間分辨率為59.5 ms，因此呼氣末期包含一個心動周期圖像數可以通過下式計算：
| ${N_p}=\frac{{Tc}}{{Tt}}$ |
其中Tc表示心動周期，Tt表示時間分辨率，Np表示一個心動周期圖像數。在每層圖像序列中，呼氣末期左心室中心相對位置高于吸氣末期，可以通過計算最大Np連續圖像相對位置的和確定呼氣末期。用N表示每層圖像數，呼吸信號可以劃分為N-Np+1區域，每個區域可以表示為：
| $\sum\limits_{t=1}^{{N_p}} {H\left(t \right), } \sum\limits_{t=2}^{{N_{p + 1}}} {H\left(t \right), } \cdots, \sum\limits_{t=N-{N_p} + 1}^N {H\left(t \right)} $ |
H(t)表示位置t之間相對位置之和。通過式(1)、式(2)獲得心動周期包含圖像數目，以及該心動周期在呼氣末期的起始位置。如圖 5所示，顯示了Np圖像從A點到B點的呼氣末期區間。
2.4 ED和ES圖像自動定位
大多數文獻中ED和ES圖像檢測依賴于左心室血池精確分割，其過程計算量大且容易受噪聲影響。本文采用文獻[4]的方法，即計算呼氣末期區間圖像之間最小互相關系數來確定ED和ES圖像。由于實時心臟磁共振成像每層圖像心臟位置不同，為了降低呼吸運動干擾，從左心室中心擴展選取一個包含整個左心室的方形區域計算互相關系數，確定ED和ES圖像，最后通過橢圓擬合面積判斷ED和ES先后順序。
2.5 統計學處理
為了驗證檢測結果是否正確，將每層檢測出的ED和ES圖像導入QMass MR(Medis，Netherland)軟件進行左心室功能分析，包括心肌質量(myocardial mass, MMASS)，舒張末期容量(end-diastole volume, EDV)，收縮末期容量(end-systole volume, ESV)，每搏血量(stroke volume, SV)和射血分數(ejection fraction, EF)參數，乳突肌不包含在心肌質量當中。對于每個參數，均計算其均值差及百分比。雙尾成對t檢驗統計差異評估，P＜0.05，認為差異有統計學意義。
3 結果
標準成像序列和實時成像序列獲得的短軸位ED和ES圖像，如圖 8所示。兩種成像序列的左心室功能參數比較如表 1所示，統計結果如表 2所示。結果表明兩種成像序列獲得的MMASS，EDV，ESV，SV，EF數據相似，差異不具統計學意義(P＞0.05)。
圖8
兩種成像序列ED和ES圖像比較
Figure8.
Comparison of ED and ES images between two imaging protocols
表1
兩種成像序列左心室功能參數比較
Table1.
Comparison of LV function analysis between two imaging protocols
表2
雙尾成對t檢驗統計結果
Table2.
Results of two-tailed paired t test
4 討論及結論
本文為了避免采用不同位置心臟圖像進行數據分析不準確的問題，通過檢測每層呼氣末期ED和ES圖像，同標準成像方法圖像進行比較。實驗結果表明應用本文方法，標準成像和實時心臟磁共振成像數據的差異不具有統計學意義，說明該法能夠正確檢測出實時心臟磁共振成像呼氣末期ED和ES圖像用于心臟功能分析。而且本法只需較少人工操作，即可自動獲得呼吸運動曲線以及ED和ES圖像。同傳統的人眼觀察挑選出ED和ES圖像方法相比，減少了后處理時間，分析人員可以直接對自動獲得的ED和ES圖像運用分析軟件進行功能分析，在下一步研究中，課題組將解決實時心臟磁共振成像圖像左心室精確分割問題。
0 引言
心臟磁共振成像是一種精準的可重復用于左心室(left ventricle, LV)功能分析的方法[1]。目前，實時自由呼吸心臟成像可用來進行左心室的功能分析。與標準的屏氣穩態自由進動序列(steady-state free precession, SSFP)相比，實時成像能縮短數據采集時間且不需要屏氣，有助于不能屏氣的患者[2]。
然而，實時自由呼吸心臟成像左心室功能分析非常耗時，因為需要從多幅磁共振序列圖像中定位心臟舒張末期(end-diastole, ED)和收縮末期(end-systole, ES)圖像。并且由于呼吸運動，左心室的位置隨著橫膈膜的運動而發生變化，如果采用每層不同心臟位置的圖像進行數據分析，可能會導致左心室功能分析數據出錯, 其主要原因是橫膈膜運動導致每層心臟處于不同位置而產生形變。為了減少功能分析誤差，選擇相同心臟位置圖像，能有效解決這個問題。而呼氣末期相對吸氣階段，橫膈膜運動較平緩，因此可以使用左心室位于呼氣末期的ED和ES圖像分析心臟功能。目前，ED和ES圖像的選取主要采用人眼觀察的方法，從每層心臟磁共振序列中找出ED和ES圖像[3]。這種方法增加了后處理時間且非常費力，因此如果可自動定位ED和ES圖像，有利于節省數據處理時間，便于后續數據分析。當前在一個心動周期自動定位ED和ES圖像的方法[4]，需要受試者屏氣，而實時心臟磁共振成像中橫膈膜的運動導致每層心臟位置不同，因此無法直接采用該方法進行測定。
本文旨在提供一種能自動檢測實時心臟磁共振成像中ED和ES圖像的方法。由于心臟的位置隨著橫膈膜的運動而發生變化，左心室的運動軌跡可以通過呼吸運動信號進行擬合, 然后運用互相關系數檢測呼氣末期區域ED和ES圖像。利用本文方法僅需較少人工干預，即可精準檢測出ED和ES圖像，并用于左心室功能分析。
1 數據來源
本研究從中國科學院深圳先進技術研究院招募10名健康志愿者，年齡為(25±4)歲，無既往心臟病史，所有受試者均簽署知情同意書，并經過單位倫理審查，同意試驗。研究使用中國科學院深圳先進技術研究院3.0 T TIM TRIO西門子磁共振儀進行掃描。首先使用標準二維SSFP的具有心電圖(electrocardiogram，ECG)觸發和呼吸抑制的序列掃描心臟，10個短軸位層包含心臟從心尖到基底；實時自由呼吸的SSFP序列在短軸位對相同層位置進行掃描，并在時間方向運用K-L變換(Karhunen-Loeve Transform，KLT)濾波增加信噪比[5]。成像參數為：層厚8 mm，間隙2 mm，視野為340×287 mm2。對于標準SSFP心臟成像序列：重復時間/回波時間=3.4/1.5 ms，矩陣大小為256×216，并行采集(integrated parallel acquisition techniques, iPAT)加速因子為2，帶寬977 Hz/像素，時間分辨率40.7 ms。對于實時心臟成像序列：重復時間/回波時間=2.5/1.1 ms，矩陣大小為160×128，時間并行采集(temporal parallel acquisition techniques，TPAT)加速因子為4，帶寬1 488 Hz/像素，時間分辨率59.5 ms，并且掃描每層持續時間為5 s，包含84幅從呼氣末期到吸氣末期心臟圖像。
2 方法
2.1 左心室定位和分割
左心室定位首先需要確定心臟在磁共振圖像上的位置。以往定位心臟的方法主要針對屏氣心臟成像，并且假設心臟是唯一運動的器官[6-8]，而實時成像中，受呼吸運動影響，心臟不是唯一的運動器官，因此前述方法不能用于實時心臟成像圖像的心臟定位。文中采用人工定位的方法，對每層所有圖像求和，手動選取求和后圖像中左心室血池中心點或者離中心距離最近較亮的點做為種子點，如圖 1所示。圖 1左圖為1名志愿者心臟短軸位求和后圖像，右圖為種子點選取圖像，該點通過手動點擊獲取，并記錄坐標。以往左心室心內輪廓分割方法采用最優閾值，能有效實現左心室血池分割[9-10]。本研究同樣采用Otsu方法[11]計算每幅圖像分割閾值，將種子點作為每幅圖像區域增長的中心，相鄰點高于閾值加入種子點，低于閾值作為心內輪廓。通過該法能快速分割每層圖像左心室血池。左心室ED和ES血池分割結果如圖 2所示。
圖1
心臟磁共振短軸位求和圖像
Figure1.
Total image of short axis of cardiac magnetic resonance
圖2
左心室血池分割結果
Figure2.
Segmentation result of LV blood pool
2.2 橢圓擬合和呼吸運動擬合
為了降低乳突肌和肌小梁對于左心室中心檢測的影響，橢圓擬合之前先對分割血池進行凸包處理, 獲得最小多邊形邊界。如圖 3所示，左圖為左心室血池分割結果，右圖為運用凸包處理獲得的邊界。然后運用最小二乘法對左心室橢圓擬合，定位左心室中心。橢圓擬合圖如圖 4所示，其中藍色為擬合的橢圓，綠色為凸包邊界，紅色為左心室中心。
圖3
血池和凸包邊界
Figure3.
Blood pool and convex hull border
圖4
ED和ES橢圓擬合
Figure4.
Ellipse fitting of ED and ES
由于心臟位置隨著橫膈膜運動而變化，左心室在呼氣末期相對位置高于吸氣末期，可用心臟空間位置描述呼吸運動軌跡。然而由于心臟跳動，信號曲線不平滑，且考慮到心臟跳動頻率遠高于呼吸頻率，因此采用低通濾波器平滑曲線以獲得呼吸運動信號，漢寧窗降低信號波動。如圖 5所示，為1名志愿者呼吸運動擬合曲線，其中橫坐標表示每層心臟圖像序列，縱坐標表示左心室中心相對圖像邊界距離，心臟和呼吸運動信號用不同符號在圖中標記，AB區間為呼氣末期區間，包含一個心動周期。
圖5
呼吸運動擬合曲線
Figure5.
Fitting curve of respiratory motion
在心尖位置，由于左心室血池面積較小，初始種子點位置固定，呼吸運動可能導致收縮末期種子點位于血池之外，如圖 6所示。為了獲得呼吸運動信號，只對種子點位于血池內圖像進行橢圓擬合，記錄中心點位置，而血池外點排除，其中心位置通過周圍點插值獲得。最后通過多次信號擬合獲得呼吸運動信號。而在基底層，當心肌包圍左心室腔體面積超過50%時，進行橢圓擬合定位左心室中心，其它圖像的左心室中心位置同樣通過插值擬合獲得。基底層連續四個心臟周期的左心室ED和ES擬合圖像如圖 7所示，其中藍色為橢圓擬合曲線、紅點為左心室中心和綠色為凸包邊界。
圖6
左心室心尖收縮末期
Figure6.
ES of LV apex slice
圖7
基底層ED和ES擬合結果
Figure7.
Fittings results of ED and ES in basal slice
2.3 確定呼氣末期
定位ED和ES圖像之前需要確定呼氣末期位置。在本研究中，心跳范圍為60~100次/min，時間分辨率為59.5 ms，因此呼氣末期包含一個心動周期圖像數可以通過下式計算：
| ${N_p}=\frac{{Tc}}{{Tt}}$ |
其中Tc表示心動周期，Tt表示時間分辨率，Np表示一個心動周期圖像數。在每層圖像序列中，呼氣末期左心室中心相對位置高于吸氣末期，可以通過計算最大Np連續圖像相對位置的和確定呼氣末期。用N表示每層圖像數，呼吸信號可以劃分為N-Np+1區域，每個區域可以表示為：
| $\sum\limits_{t=1}^{{N_p}} {H\left(t \right), } \sum\limits_{t=2}^{{N_{p + 1}}} {H\left(t \right), } \cdots, \sum\limits_{t=N-{N_p} + 1}^N {H\left(t \right)} $ |
H(t)表示位置t之間相對位置之和。通過式(1)、式(2)獲得心動周期包含圖像數目，以及該心動周期在呼氣末期的起始位置。如圖 5所示，顯示了Np圖像從A點到B點的呼氣末期區間。
2.4 ED和ES圖像自動定位
大多數文獻中ED和ES圖像檢測依賴于左心室血池精確分割，其過程計算量大且容易受噪聲影響。本文采用文獻[4]的方法，即計算呼氣末期區間圖像之間最小互相關系數來確定ED和ES圖像。由于實時心臟磁共振成像每層圖像心臟位置不同，為了降低呼吸運動干擾，從左心室中心擴展選取一個包含整個左心室的方形區域計算互相關系數，確定ED和ES圖像，最后通過橢圓擬合面積判斷ED和ES先后順序。
2.5 統計學處理
為了驗證檢測結果是否正確，將每層檢測出的ED和ES圖像導入QMass MR(Medis，Netherland)軟件進行左心室功能分析，包括心肌質量(myocardial mass, MMASS)，舒張末期容量(end-diastole volume, EDV)，收縮末期容量(end-systole volume, ESV)，每搏血量(stroke volume, SV)和射血分數(ejection fraction, EF)參數，乳突肌不包含在心肌質量當中。對于每個參數，均計算其均值差及百分比。雙尾成對t檢驗統計差異評估，P＜0.05，認為差異有統計學意義。
3 結果
標準成像序列和實時成像序列獲得的短軸位ED和ES圖像，如圖 8所示。兩種成像序列的左心室功能參數比較如表 1所示，統計結果如表 2所示。結果表明兩種成像序列獲得的MMASS，EDV，ESV，SV，EF數據相似，差異不具統計學意義(P＞0.05)。
圖8
兩種成像序列ED和ES圖像比較
Figure8.
Comparison of ED and ES images between two imaging protocols
表1
兩種成像序列左心室功能參數比較
Table1.
Comparison of LV function analysis between two imaging protocols
表2
雙尾成對t檢驗統計結果
Table2.
Results of two-tailed paired t test
4 討論及結論
本文為了避免采用不同位置心臟圖像進行數據分析不準確的問題，通過檢測每層呼氣末期ED和ES圖像，同標準成像方法圖像進行比較。實驗結果表明應用本文方法，標準成像和實時心臟磁共振成像數據的差異不具有統計學意義，說明該法能夠正確檢測出實時心臟磁共振成像呼氣末期ED和ES圖像用于心臟功能分析。而且本法只需較少人工操作，即可自動獲得呼吸運動曲線以及ED和ES圖像。同傳統的人眼觀察挑選出ED和ES圖像方法相比，減少了后處理時間，分析人員可以直接對自動獲得的ED和ES圖像運用分析軟件進行功能分析，在下一步研究中，課題組將解決實時心臟磁共振成像圖像左心室精確分割問題。
