1.寫在前面
這個項目從我個人的SRT（student research training）演變而來，研究目的是對豬心左心室磁共振PSIR圖像心肌梗死區域進行分割，由於水平和數據有限，該項目不涉及機器學習等方法，而着重於一些傳統方法，希望對大家有所幫助。
github項目網址https://github.com/BertramRay/CMR_postprocess_training
github上下載完所需數據和代碼，簡單配置后即可嘗試快速運行。
2.摘要
在本項目中，我們從圖像分割閾值，黑核識別，噪聲去除等多個不同的角度入手，探究並初步獲得了一個利用傳統方法對心臟磁共振圖像心肌梗死區域進行分割的整合方案。我們小組1.討論並嘗試了n-SD法，FWHM法，雙峰法等圖像分割閾值確定方法。2.引入了生長法去除由於血池導致的ROI邊緣區域異常高亮信號。3.引入了種子點對圖像梗死區域范圍進行限制，以消除由於磁場不均勻性導致的遠端高亮偽影。4.提出並實現了黑核閉合（半閉合）狀態檢測方法。5.加入了部分距離加權的梯度信息以提高黑核檢出率。最后，我們小組將上述方法進行整合，並在多個不同病例上驗證了該方案的普適性。
3.介紹
3.1. 相位敏感反轉恢復序列Phase sensitive inversion recovery (PSIR)
圖 1 IR和PSIR圖像亮度與Mz的關系1
相較於IR序列得到的圖像，PSIR序列得到的圖像能夠區分出不同極性相同大小的磁化矢量，從而在某些情況下提供更加豐富的對比度信息。
圖 2 典型的PSIR序列2
為了區分不同極性的磁化矢量，需要引入參考采樣脈沖。
3.2. 釓劑延遲增強Late gadolinium enhancement (LGE)
利用LGE區分正常心肌與梗死心肌利用了釓元素能夠顯著縮短T1值這一特性，相比於正常心肌，梗死心肌的T1加權圖像的信號尖峰出現較晚，所以在梗死心肌信號最大值處采集所得的圖像表現為梗死區域高亮而正常心肌偏暗。
圖 3 LGE區分不同心肌的原理示意圖3
4.具體工作
4.1. 圖像分割閾值的選定
4.1.1.n-SD法
通過手動選定種子點，確定人體正常肌肉組織參考區域，計算參考區域內像素點灰度值的平均值與標准差，有
我們認為，灰度值高於該閾值的像素點為梗死心肌組織，而灰度值低於該閾值的像素點為正常心肌組織或黑核。一般來說，n設置為5。該方法在所有方法中表現最為穩定。
圖 4 n-SD法分割得到的梗死區域（右上）FWHM法分割得到的梗死區域（右下），可見FWHM法所獲得的閾值偏低
4.1.2.FWHM法
在此基礎上我們還考慮了FWHM法，即首先統計參考心肌組織區域各灰度值像素點數量的分布情況，尖峰右側半峰高所對應的灰度值即為所求閾值。該方法偶見所得閾值偏低的情況。
4.1.3.雙峰法
該方法假設了梗死區域和正常心肌組織區域的像素點統計分布呈雙峰，並尋找雙峰間的低谷，低谷對應的灰度值即為所求的閾值。該方法在實際測試中表現不夠理想，原因在於高亮梗死區域像素點灰度值分布並不呈現尖峰。如圖所示，縱坐標為像素點個數，橫坐標為灰度值，高亮梗死區域像素點灰度值跨度較大且未見明顯尖峰。
圖 5 ROI區域像素點灰度值分布情況
事實上，絕大多數基於雙峰假設的圖像閾值確定方法在該分割任務中表現均不理想，如OTSU法，迭代法等。因此最終我們決定在后續梗死區域圖像分割任務中統一采用n-SD方法。
4.2. 基於生長法的血池導致的邊緣異常高亮的去除
在實際使用n-SD法進行梗死區域圖像分割時，我們發現由於前期ROI區域較大，導致部分血池被划入到ROI區域中，而血池的信號強度較大，易與梗死區域混淆。為了降低血池帶來的干擾，我們利用血池邊緣高亮像素點較少且不連續的特點，設計並編寫了生長法的代碼，該方法能夠識別計算圖像中各個高亮區域的高亮像素點個數，低於某一特定值的高亮區域將被填充為黑色背景區域，而成片連結的梗死區域則不受影響。一般來說，該任務中該像素點個數閾值設定為6-12取得的效果最為理想，如下圖所示。
圖 6 不同強度的生長法去除邊緣異常高亮的效果
4.3.閉合黑核與半閉合黑核的識別與修正
黑核是一類特殊的心肌梗死區域，相比於正常的高亮梗死區域，黑核有如下特點：1.灰度值顯著低於高亮梗死區域，與正常心肌灰度值大致相當或更低。2.大多數情況下黑核被高亮梗死區域包圍，少部分情況下黑核緊貼心內膜存在。3.包圍黑核的高亮梗死區域在圖像上可能是閉合的，也可能是近似閉合的，視閾值與梯度加權信息而定。
為識別並填充黑核，我們首先編寫了能夠識別閉合黑核的代碼，該方法能夠利用連通性，將背景黑色區域與黑核黑色區域區分開來，並將黑核黑色區域填充為梗死區域。如下圖所示，左圖中的黑核被成功填充。
圖 7 閉合黑核的識別示例圖
為填充未被閉合的黑核，我們提出了一種邊緣試探的方法。如下圖所示，我們將所有邊緣白色高亮像素標記為深藍色方塊，並令其自然延伸出一個邊緣試探的淺藍色方塊，再令淺藍色方塊延伸出三個綠色的基礎檢測區域方塊，若該基礎檢測區域與淺藍色或深藍色方塊重疊，則視為試探閉合成功，此時將淺藍色方塊標記為高亮梗死區域。同時，還可以設定一定長度的擴充檢測區域，當橙黃色的擴充檢測區域與邊緣試探像素重合時，視為試探閉合成功，填充所有橙黃色方塊。隨后再次進行閉合黑核填充，此時原本未閉合的黑核也會變成閉合的黑核從而被識別。效果如下圖所示。
圖 8 邊緣試探法填充未閉合黑核示例圖
注意到擴充檢測區域的橙黃色方塊數量直接決定了能夠閉合的黑核的數量，但是過多的擴充檢測區域可能會使高亮梗死區域的邊緣信息部分喪失，在實際使用時還需注意平衡。
4.4.利用種子點消除遠端磁場不均勻偽影
在實際心臟磁共振操作中由於操作不當或是受線圈和儀器的影響，會出現磁場不均勻的情況，進而導致最終圖像大范圍區域內的灰度值異常升高或降低。為盡可能減少磁場不均勻性導致的偽影，我們引入了梗死區域種子點。
種子點能夠保留與之相接或范圍內的高亮梗死區域，使得出現在遠端的由於磁場不均勻導致的高亮偽影能夠被填充為黑色背景區域。
圖 9 遠端磁場不均勻偽影移除示例圖
4.5.考慮梯度信息對圖像做銳化處理
通過對ROI區域中高亮梗死區域，黑核區域和正常心肌區域的灰度值像素點進行統計后，我們發現，正常心肌與黑核灰度值有較大重合區域，難以通過閾值方法簡單划分，而且黑核的識別檢測嚴重依賴周圍的梗死區域，如果黑核周圍的梗死區域灰度值低於一般的高亮梗死區域，那么黑核的檢出就會受到嚴重的影響。
如下圖所示，右側圖像紅色部分為應檢測出的高亮梗死區域，而藍色部分為黑核區域與灰度值偏低的高亮梗死區域，在考慮梯度信息前難以識別出黑核和較暗的高亮梗死。
圖 10 實際黑核情形示例圖
為此，我們在圖像處理的最開始引入梯度信息並與原圖像信息進行結合。即計算ROI內每個像素點與周圍四個像素點的灰度值差的絕對值的平均值，並於原像素點灰度值進行加權。這樣原圖像的灰度值陡變區域表現為更高的灰度值，能夠有效的提高黑核被閉合的程度。具體結果如下圖所示。
圖 11 不考慮梯度信息（左）考慮灰度信息（右）
5. 數據計算
5.1. PSIR多層面心肌梗死比例計算
通過分別計算不同層面的梗死區域和ROI區域面積，再累加相比，能夠粗略得到實驗對象的心肌梗死比例。
| 實驗對象 |
TotalVolumn /pixel |
LGEVolumn /pixel |
LGERatio |
| C139 |
66027 |
17264 |
26.15% |
| B306 |
40463 |
9018 |
22.29% |
表1 不同實驗對象的心肌梗死比例計算
經驗證，該比例與實際梗死區域比例大致相當。
5.2. 不同病例的方法普適性驗證
經驗證，該方法在不含黑核或黑核較少的實驗對象中結果穩定性好，程序運行速度快。在含黑核較多的實驗對象中，普遍存在靠心內膜的黑核難以檢出的問題。
| 實驗對象 |
含梗死區域的全部圖像數 |
完全檢出黑核的圖像數 |
黑核完全檢出率 |
| C139 |
29 |
25 |
86.21% |
表2 C139實驗對象黑核檢出率計算表
注：B306心肌梗死區域內無黑核
6. 項目總結與致謝
6.1.PSIR圖像梗死區域分割程序處理流程總結：
一.數據導入與參數設置。
二.種子點選擇（梗死中心區域，n-SD參考區域）
三.梯度信息處理
四.n-SD法處理
五.邊緣噪點處理
六.閉合法與半閉合法處理黑核
七.遠端磁場偽影處理
八.多層面圖像梗死比例計算
6.2.致謝
參考文獻
[1] Question and Answer in MRI - Magnitude VS Phase Sensitive IR Allen D. Elster, MD FACR
http://mriquestions.com/phase-sensitive-ir.html
[2] Kellman P et al. Magn Reson Med. 2002;47:372–83.
[3] Kellman P, Arai AE. J Magn Reson Imaging. 2012;36:529–42.