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.