卷積神經網絡結合可變模型的MRI圖像自動分割方法

任 俠 胡玉平

1(衢州職業技術學院信息工程學院 浙江 衢州 324000)2(廣東財經大學信息學院 廣東 廣州 510320)

0 引 言

磁共振成像(Magnetic Resonance Imaging，MRI)技術是一種在現代醫學領域中應用最為廣泛的成像技術[1-2]。在臨床中，基于心臟MRI數據庫的左心室分割有助于精確算心室容積、射血分數、左心室質量和壁厚等關鍵指標，以及分析壁厚異常等情況[3-4]。傳統的心室分割主要依賴于醫學專家手繪與臨床實踐，不僅浪費時間，而且容易發生偏差。因此需要開發左心室自動分割方法，以便加速診斷進程。在此過程中，圖像分割[5-6]占據著非常重要的位置，其性能的優劣直接關系著診治的效果。左心室MRI自動分割存在著左心腔異構、心肌信號強弱變化、噪音等多個難題，使得實現基于MRI的心室自動分割成為一個極具挑戰的研究課題。

心室自動分割方法主要包括：像素分類法、圖像法、形變法、圖集模型法。像素分類法、圖像法和形變法魯棒性和精確度較低，并且需要借助廣泛的用戶交互。圖集模型法的魯棒性和精度較高，并且通過建立通用模型可以明顯減少用戶交互。然而，如何建立能夠完全包含心臟所有可能的形狀以及動態變化的模型，是一個異常復雜的問題。如果選用的數據庫較小，就會導致心臟分割產生較大偏差，當心臟的形狀超出了學習的范圍時，會使得這種方法的效率大大降低。此外，現有的心室分割學習方法也具有一定的局限性，如隨機植入方法等通過圖像強度將分割問題定義為分類任務。這些方法利用多個階段進行強度標準化、評估和正規化，導致算法的響應速度非常慢，并會影響后續步驟的成功率。

針對上述問題，文獻[7]利用各向異性擴散方法對心臟核磁共振圖像進行濾波處理，進而利用簡化脈沖耦合深度學習提取心內膜信息，采用形態學處理器獲取心外膜信息。文獻[8]提出了基于深度全卷積深度學習架構，同時借助端到端的模型訓練來實現MRI圖像的像素級分類。文獻[9]開發了一種LV預測方法，借助深度學習技術和大規模心臟MRI數據集實現了LV容量的準確預測，為進一步的臨床應用奠定了基礎。文獻[10]設計了兩個卷積深度學習，一個用于LV定位，另一個用于確定心臟半徑，在100個數聚集上的實驗證明了所提方法的有效性。

受上述方法的啟發，同時考慮到MRI在臨床中需要專家進行手動操作的事實，本文提出一種基于卷積神經網絡和可變模型算法的左心室MRI全自動分割方法，獲得了較高的分割精度。其主要創新點為：

(1) 在左心室位置檢測中，通過引入稀疏自動編碼技術實現卷積神經網絡的預訓練，在提高數據可用性的基礎上，避免了過擬合，實現了位置的高效精確檢測。

(2) 在左心室形狀推斷中，通過搭建并訓練堆棧稀疏編碼器來獲取左心室的形狀組合，在提高分類精度的基礎上簡化了計算復雜度。

(3) 在左心室圖像分割中，通過最小化能量方程將可變模型變為動態輪廓，在阻止輪廓的收縮或泄露的前提下，實現了圖像的精確分割。

實驗結果表明，提出的方法可以大幅度提高左心室MRI圖像分割準確度和魯棒性，實現圖像的全自動精確分割。

1 分割方法設計

圖1為本文方法的步驟。首先，利用卷積網絡對輸入圖像進行訓練，確定LV的位置以及包含LV的區域；其次，利用堆棧自動編碼器推斷LV的形狀，描繪LV的輪廓；最后，對推斷出的形狀進行初始化，并將其并入可變模型，實現圖像分割。其中，LV的輪廓定位可以減少三維重建中切片之間的偏差，每個步驟分別進行離線訓練以獲得最優參數值，訓練后自動執行圖像分割。

圖1 分割方法步驟

1.1 自動檢測

原始心臟MRI圖像數據庫中包括心臟及其周邊器官的信息[11]，為了降低計算的復雜度，提高分割精度，需要首先確定LV及其領域的位置。圖2為基于卷積網絡的LV位置自動檢測流程。

圖2 MRI數據庫LV自動檢測流程圖

(1) 通過將原始圖像尺寸由64×64降低至32×32并作為算法輸入，可以有效降低計算復雜度。

(2) 將過濾器Fl∈R8×8，b0∈R100與輸入圖像進行卷積運算，得到卷積特征圖像。

定義灰度值圖像I：Ω→R，Ω?R2，其像素大小為64×64，坐標為[i,j]，I[i,j]為像素強度。卷積特性的計算公式為：

(1)

式中：ζ∈[0,1]為特征系數。Zl[i,j]的計算公式如下：

(2)

式中：1≤i,j≤28；l=1,2，…,50；x[i,j]表示矩陣x的第i行、第j列；x[i]表示向量x的第i個元素。

(3) 利用平均采樣方法處理卷積特征圖像，此時可以得到3×3非重疊區域內每個特征圖像的平均值，其計算公式為：

(3)

式中：1≤i1,j1≤3；ε∈[1, 2]為校正參數。

(4) 將提取的混合特性展開為向量p∈R900，通過連接具有256個輸出信號的邏輯回歸層以產生大小16×16的LV領域。

經過以上步驟，可將大小為64×64的原始圖像降維成大小為16×16的算法輸入圖像，從而完成了圖像數據的超抽樣，進而通過中心計算將原始圖像縮減為50×50大小的鄰域。

最關鍵的為通過訓練卷積網絡獲得過濾器Fl和b0等關鍵參數。通過稀疏自動編碼獲得過濾器，進而隨機初始化過濾器，并進行大量的卷積網絡訓練，即可得到滿意的結果。稀疏自動編碼充當預訓練的角色，在有限的數據下對網絡進行訓練，同時可以避免過度擬合。

圖3展示了稀疏編碼器的訓練過程。定義W1∈R50×8為隱藏層與輸出層的權值，W2∈R8×50為輸入層和隱藏層之間的權值。隱藏層、輸出層的計算公式分別為：

a(i)=f(W2·x(i)+b2)

(4)

y(i)=f(W3·a(i)+b3)

(5)

式中：f(x)=1/(1+e-x)是S型激活函數；b2、b3為偏置向量；W3為校正權值。

圖3 稀疏編碼器訓練過程

首先，稀疏編碼器通過最小化成本函數來獲取最優解：

(6)

其次，通過最小化式(7)預訓練輸出層：

(7)

式中：U(i)為標記數據；N2為訓練數據的數量。

最后，通過最小化式(8)調整網絡：

(8)

1.2 形狀推斷

圖4為基于訓練堆棧稀疏編碼器的MRI圖像形狀推斷的原理圖。堆棧稀疏編碼器包括一個輸入層、兩個隱藏層和一個輸出層[11-12]。首先，將1.1節中得到的子圖像進行采樣、展開為向量xs∈R2 048，并將其作為輸入量送入輸入層；其次，在隱藏層通過計算h1=f(W6xs+b4)和h2=f(W5h1+b5)構建抽象特征；最后，在輸出層計算ys=f(W6h2+b6)來產生二進制圖像，且這些二進制圖像在LV邊界之外都要歸零；W4∈R50×2 048、b4∈R50、W5∈R50×50、b5∈R50、W6∈R2 048×50、b6∈R2048均為可訓練矩陣和向量。

圖4 MRI圖像形狀推斷原理

堆棧稀疏編碼器的訓練具體分為兩個步驟：預訓練和調整。由于在應用中標記數據的存取受到限制，因此需要采用離散層預訓練方法。離散層預訓練可以阻止過度擬合，通過無標記數據層層獲得堆棧稀疏編碼器的參數W4、W5，而W6可以通過標記數據獲得。具體方法如下：

(1) 從堆棧稀疏編碼器中分離出輸入層和隱藏層H1，通過將與輸入層相同尺寸的輸出層添加到輸入層和H1，即可建成稀疏編碼器，進而在無監督模式下對稀疏編碼器訓練得到W4。稀疏編碼器的輸入/輸出數據是大小為50×50的子圖像，位于LV中心，從標準尺寸的訓練圖像中提取得到。輸入圖像大小為32×32，與堆棧稀疏編碼器的輸入尺寸相匹配。完成第一個稀疏編碼器的訓練后，其輸出層即被廢棄。稀疏編碼器中隱藏單元的輸出層用作下一個隱藏層H2的輸入。

(2) 隱藏層H1和H2從堆棧稀疏編碼器中分開，構建第二個堆棧稀疏編碼器。對第二稀疏編碼器進行訓練得到W5，并且此時不需要標記數據。

(3) 將隱藏層輸出作為最后一層的輸入，在監督模式下對最后一層訓練得到W6，通過成本函數來訓練最后一層：

(9)

式中：L(i)∈R2 048為標記數據，其為二進制圖像，由手動分割生成。

通過逐層預訓練，為生成的參數W4、W5、W6選擇合適的初始值，并通過最小化式(10)的成本函數對搭建的結構進行調整。

(10)

1.3 分割和定位

利用可變模型和推斷出的形狀組合對心臟精確分割。通過最小化能量方程將可變模型變為動態的輪廓，當輪廓位于實際物體的邊緣時，能量方程達到最小值。在大多數傳統的可變模型中，由于左心室乳頭肌的存在及其與周圍組織對比不明顯，輸出輪廓往往向內收縮或者向外泄露。通過1.1節、1.2節的前期步驟推斷出心臟形狀，再將其作為計算的初始值，從而解決這些問題。此外，通過將推斷出的心臟形狀納入到能量方程中，也可以阻止輪廓的收縮或泄露。

定義輸入的子圖像作為Is:Ωs→R，φ(x,y)為基準水平，其中(x,y)為圖像像素坐標。如果像素在輪廓之內，則返回負值，反之則返回正值。此外，定義推斷形狀的基礎水平方程為φshape(x,y)，能量方程為：

E(φ)=τ1Elen(φ)+τ2Ereg(φ)+τ3Eshape(φ)

(11)

式中：τ1、τ2、τ3為能量參數增益；Elen(φ)為長度能量方程；Ereg(φ)為區域能量項；Eshape(φ)為形狀能量項。

長度能量方程的表達式為：

(12)

區域能量項的表達式為：

(13)

形狀能量項的表達式為：

(14)

式中：δ(φ)、H(φ)和▽(·)分別為δ函數、Heaviside階躍函數和梯度運算函數；c1和c2分別是心臟輪廓外部和內部輸入圖像Is的平均值；τ1、τ2、τ3均為訓練中根據經驗得出的能力參數，分別取值為τ1=1.5、τ2=1、τ3=0.5。將可變模型組合方法得出的結果與已知訓練數據庫中相應的圖像和標簽進行對比，并對參數τ1、τ2、τ3進行調整，以獲得最好的評估指標。

可變模型建立的目標是搜尋位于實際物體邊緣的特定輪廓(表示為C*或者φ*)，該輪廓可以通過最小化能量函數來獲得:

(15)

通過式(16)的梯度下降算法求解式(15)，通過將φ設置為時間的函數，同時結合歐拉-拉格朗日方程，可得：

(16)

式中：Div(·)為散度因子。

在使用梯度下降算法時，首先需要通過推斷形狀得到初始化φ(0)，進而更新迭代公式：

(17)

式中：γ1、γ2為迭代步長。

通過對式(17)進行迭代運算可以獲得最終的輪廓C*或者φ*。迭代停止的判別標準是檢查所得解是否穩定，或者通過比較當前和之前迭代中輪廓不同長度來確定。

當對心臟重新構建三維模型時，需要重新考慮圖像之間可能存在的偏差。核磁共振成像掃描時出現的偏差主要是由于掃描時病人的呼吸和移動造成的，如果忽略這些偏差就會在心臟重建時產出鋸齒狀的不連續表面，從而產生一定的干擾。為了解決這個問題，引入二次多項式對偏差進行估計和糾正。

(18)

(19)

通過求解二次多項式獲得xi、yi：

xi=a1i2+b1i+c1

(20)

yi=a2i2+b2i+c2

(21)

式中：a1、a2、b1、b2、c1、c2是未知參數，通過求解最小均方誤差來得到:

(22)

(23)

在得到這些未知參數的估值后，使用基于線性插值的仿射變換獲得心臟輪廓，最后根據測定出來的中心坐標值即可獲得線性的輪廓堆棧。

2 實驗分析

2.1 實驗概述

為了驗證所提方法的有效性，在某AI挑戰賽的訓練數據組中收集所有情況的圖像和輪廓，該數據庫包含30個核磁共振成像(MRI)數據庫，分成三個數據庫，每個數據庫包含10個案例，這些案例包括4個缺血性心力衰竭(SC-HF-I)、2個非缺血性心力衰竭(SC-HF-NI)、2個左心室肥大(SC-HYP)和2個正常情況(SC-N)案例[12-13]。此外，在心臟舒張末端和收縮末端中，專家對心室的手動分隔也包括在數據庫中。一個典型的數據庫可以從底部到尖端每4～8個短軸切片中包含15幀，圖像參數：厚度為5毫米，圖像大小為128×128像素。利用該數據庫對本文方法進行校驗和在線評估。

在進行實驗時，需要將該數據庫分成大輪廓組和小輪廓組底部或中間附近的圖像切片屬于大輪廓，中心頂端屬于小輪廓，因為靠近頂端的輪廓要比底部的輪廓小很多，每一組分別有大約100和80幅圖像。為了擴大數據庫以便于算法的驗證，需要使用圖像轉移、圖像旋轉、更改像素強度等辦法人工擴大訓練數據組，這些技術使得訓練數據組增大十倍，最終每一個組中分別有1 150和950幅圖像或標簽。

由于需要學習大量參數，因此深度學習網絡可能會出現一定程度的過度擬合現象，為了解決或防止過度擬合的問題，需要采用層預訓練、l2規則化以及稀疏約束。雖然訓練數據缺失是一個較為嚴重的問題，但是可以通過預訓練的方法進行解決，該方法可以保持隱藏層較小的單元數，而且不超過三層，進而保證參數的數量易于掌控。在訓練過程中，通過將訓練數據集分成9個訓練主題和3個驗證主題，并借助早期停止監控，可以阻止過度擬合，最后就可以進行交叉驗證。此外，需要將訓練數據集人為放大。網絡的超參數，包括層數和單元數、過濾器數、過濾池大小等，在交叉驗證過程中可以通過經驗來確定。

至于軟硬件實驗平臺，本文選用MATLAB 2013a和聯想高精密工作站，該工作站的硬件配置為Intel(R) Xeon(R) CPU 2.6 GHz，64 GB RAM，操作系統為64位 Windows 10[14-15]。

為了表明所提算法的優越性，在精度、平均垂直距離、Dice度量、豪斯多夫距離、輪廓比例和一致性等指標上，將提出的方法與專家手動標注方法進行比較。一般認為，如果平均垂直距離小于5 mm，則表明分割方法比較好。平均垂直距離是指自動分割輪廓與專家手動注釋輪廓之間的距離，取所有輪廓點的平均值，如果該值比較高，則代表兩個輪廓匹配度較低。輪廓重疊度量為Dice度量，其取值(介于0和1之間)為DM=2Aam/(Am+Aa)，Aam為交集，Am為輪廓面積手動分割，Aa為輪廓面積自動分割。Dice度量取值越高，表明手動分割和自動分割之間的匹配度越好。豪斯多夫距離是指手動和自動輪廓之間最大的垂直距離。一致性系數是指未分割像素與正確分割像素數量之間的比例，定義為CC=(3DM-2)/DM。

此外，在自動和手動LV分割結果的基礎上，計算三個臨床參數，包括舒張末期容積(EDV)、收縮末期容積(ESV)以及射血分數(EF)[16]，并使用這三個指標來進行相關性和Bland-Altman分析。為了評估觀測器之間及其內部的變化性，需要計算變異系數(即自動和手動結果之間差異標準偏差除以平均值)。

2.2 實驗結果及分析

圖5展示了一個典型心臟MRI數據庫的手動及自動LV分割結果，以及重建LV室的三維圖(前部、底部、頂部)。

由于存在乳頭肌和低分辨率，導致LV頂部及中間的圖像切片分割結果比較復雜。圖6展示了存在擾動時，三個典型心臟MRI數據庫的手動及自動LV分割結果，自動分割結果用實線圓圈表示，專家手動分割結果用虛線圓圈表示。

圖5 心臟MRI自動與手動LV分割結果

圖6 存在擾動時心臟MRI自動與手動LV分割結果

圖7的每一行對應一個病人/數據集，(a)、(b)、(c)分別為訓練、在線和驗證三個階段，在每個子圖中由上至下分別為缺血性心力衰竭、非缺血性心力衰竭、LV肥大和正常四種心臟狀況。

(a)

(b)

(c)圖7 多種情況下自動與手動LV分割

表1為驗證數據集和在線數據集計算指標的平均值。每一個數據集的兩行分別對應推斷圖像中的原始輪廓以及最終輪廓。表2為本文方法與其他方法在相同數據集中進行圖像分割時的性能比較結果。

表1 驗證算法的評價指標

表2 提出的方法與其他方法的分割性能比較

由圖5可知，在對LV從底到上進行分割時，對齊過程產生了一個光滑的三維 LV重建圖，圖5的第一幅圖像與實地狀況相比有一點點輕微的泄露，由于邊緣的模糊性以及在可變模型中輪廓傾向于向周圍泄露，因此這種情況也是分割中的一個關鍵問題。通過將推理圖形與可變模型相結合，可以避免該圖中的泄露情況，例如圖7(c)第一行的第一幅及第二幅圖像就避免了這一情況。圖6所示的可變模型，由于乳頭肌的存在傾向于內向收縮，或由于分辨率較低及圖像頂端較小的對比度向外泄露，但是提出的方法均可以克服上述缺點。

表1中的指標說明推斷圖形提供了較好的原始輪廓，且準確度達98.2%(DM角度)，整合后的可變模型提供了最終的輪廓，其他指標也得到了明顯的提高。由表2可見，本文方法在良好輪廓比例和一致性兩個指標上優于其他算法，在Dice度量和平均垂直距離兩個指標上與現有方法相差并不大，這是因為本文方法關鍵在于以最小的復雜性實現較好的分割精度，并且盡可能地貼近專家手動分割效果。因此，本文方法在分割精度上可以滿足上述要求。

圖8-圖10展示了三個臨床心臟指數之間的相關性。手動和自動之間的高相關性表明提出的方法對于LV功能評估的準確性及臨床應用性具有重大提升。

圖8 舒張末期容積(ESV)的相關性分析圖

圖9 收縮末期容積的相關性分析圖

圖10 射血分數的相關性分析圖

通過測量運行時間可以看出，本文方法可以在一個合理的時間范圍內完成訓練任務(該過程也可離線完成)。卷積網絡所需的時間最長，并且需要包含圖像過濾器的卷積。然而，可以將算法開發成在CUP平臺進行計算來縮短時間。在測試過程中，常規圖像LV分割的平均執行時間少于0.4 s，其中大部分時間用于卷積網絡及可變模型。由于推斷圖形的初始化及整合，整合后的可變模型的收斂速度要快于傳統的可變模型。

盡管三維方法在很多醫學圖像分析中已經非常先進，但是本文提出的方法仍然采用二維處理方法。這是因為心臟MRI固有的三個問題阻止了三維分析的進一步發展：(1) 大多數通過常規方法獲取的切片之間的差異(垂直維)相對較大，切片之間的像素強度很難準確評估。(2) 由于心臟MRI中的移動偽影也導致了切片之間難以對準，從而導致不同切片的腔中心不能在同一個位置。臨床實踐中心室分割主要用于計算臨床指標，而本文方法可以提供這些參數，而且二維過程進行相應的計算精確度更高。(3) 基于深度學習的心臟MRI分割的主要缺點是缺少足夠的訓練及驗證數據，對于深度學習網絡而言，更多數據意味著更好的普適性，并可削弱過度擬合問題。然而目前并沒有相應的分析方法來設計深度學習網絡中的超參數，一般只能依賴經驗值。

3 結 語

本文提出了一種左心室自動分割方法，通過將分割問題分解為基于卷積神經網絡的定位問題、基于堆棧稀疏編碼器的圖形推理問題以及基于可變模型的圖像分割，提高了左心室圖像分割的邏輯性與精度。實驗結果表明，本文方法具有較高的分割精度。

未來將進一步研究基于深度學習的左心室圖像分割方法，借助深度學習方法超強的特征提取能力來進一步優化定位、推理與分割問題。