基于FCM聚类算法与自适应阈值的冠脉钙化分割与量化

赵 聪 陈晓冬∗ 张佳琛 汪 毅 贾忠伟 陈向志

1（天津大学精密仪器与光电子工程学院光电信息技术教育部重点实验室，天津 300072）

2（中国人民解放军第254医院心血管内科，天津 300142）

3（中国人民解放军第254医院放射科，天津 300142）

基于FCM聚类算法与自适应阈值的冠脉钙化分割与量化

赵 聪1陈晓冬1∗张佳琛1汪 毅1贾忠伟2陈向志3

1（天津大学精密仪器与光电子工程学院光电信息技术教育部重点实验室，天津 300072）

2（中国人民解放军第254医院心血管内科，天津 300142）

3（中国人民解放军第254医院放射科，天津 300142）

针对基于CTA图像进行冠脉钙化量化时存在的无法克服噪声以及阈值选择不稳定问题，提出一种基于聚类算法与自适应阈值的冠脉钙化分割与量化方法。首先根据CT值和空间位置对冠脉血管内的像素点构建特征向量，继而根据血管骨架点数目构建自适应聚类数，使用模糊C均值（FCM）聚类算法将冠脉区域划分为CT值分布相似的区域；然后使用高斯函数拟合冠脉灰度直方图，根据高斯拟合参数构造自适应阈值，对上述区域进行钙化分割；最后根据分割结果，参考Agatston钙化分量化标准进行钙化分计算。在30组人体冠脉CTA数据的测试结果中，对冠脉钙化量化的灵敏度和特异性分别达到89.5%与98.6%，计算得到的钙化体积和Agatston钙化分与标准结果的皮尔逊系数分别为0.974与0.975，远高于同类型基于一阶微分进行阈值选择方法（DBTD）对应的0.523与0.501。实验结果表明，该方法可用于冠脉钙化分割与量化，且具有全自动、鲁棒性好、能有效抗噪等特点。

CT图像处理；冠脉钙化；分割与量化；模糊C均值（FCM）聚类；钙化分计算

引言

冠心病是指冠状动脉血管由于发生粥样硬化、炎症、栓塞等病变而引起血管腔狭窄或阻塞，进而造成心肌缺血、缺氧或坏死的现象，特别是当动脉硬化斑块突发破裂时可能引起心血管突发性堵塞，进而造成病人猝死，是世界上致死率最高的疾病之一，严重威胁着我国国民的生命健康安全［1-3］。医学研究表明，对冠脉钙化的探测和量化会为预测冠心病的发展提供重要的风险评估依据［4-6］，从 CT（computed tomography）数据中获得的冠脉钙化分已经被证明是冠心病的有效预测参数［6］，而钙化斑块的形状、大小与位置也是医生在采取介入式治疗手段时需要参考的重要信息［7］。近年来，越来越多的团队开始研究在 CT造影（computed tomography angioplasty，CTA）数据中分割、量化冠状动脉钙化并计算钙化分的方法，这不仅能够减少患者接受扫描次数和受到辐射的剂量［8］，而且由于CTA的分辨率一般高于CT数据，因此也可以提高钙化量化精度。一些研究工作通过机器学习等方法对比CT与CTA数据中的钙化分，并获取两者之间的经验转化系数，进而得到在 CTA中量化钙化的方法［9-11］，但这些方法的鲁棒性往往受到CT协议、CT扫描仪生产厂商、造影剂用量等众多因素的影响［11］。

为了克服上述问题，科研工作者们提出了一些自适应的算法。Eilot等提出一种基于模型的钙化量化方法［12］，通过虚拟移除造影剂进行钙化探测和钙化分计算，虽然对小型钙化斑块有着较好的识别能力，但执行过程复杂、运算消耗大。Wolterink等提出了一种基于模式识别的钙化识别方法［13］，利用神经网络和随机森林算法进行钙化探测，其训练过程中需要大量样本，且对硬件要求高。Teßmann等提出了一种基于直方图一阶微分阈值选择（derivative-based threshold determination，DBTD）的钙化分割和量化方法［14］，得到的阈值不仅与人工选择的阈值相对接近，且具有相对较快的计算速度。然而，该方法依赖于相对光滑的直方图和经验参数，当CTA数据存在较多噪声或直方图波动强烈时，其鲁棒性将不可避免地受到影响；而且由于造影剂的使用，靠近冠脉开口处的图像噪声的CT值可能已经超过了局部钙化斑块，因此单纯依靠全局阈值对钙化进行探测可能会导致误分割或欠分割问题。

本研究针对钙化探测及量化时存在的无法克服噪声与阈值选择不稳定问题，提出了一种基于聚类算法与自适应阈值的冠脉钙化分割与量化方法。该方法先通过模糊C均值（fuzzy C-means，FCM）聚类算法对冠脉血管进行区域划分，继而使用高斯函数拟合冠脉灰度直方图以构建自适应阈值进行钙化分割，最后根据筛选结果计算钙化分。由于血管区域划分和自适应阈值的使用，该方法具有较高的灵敏度和特异性，且能够有效地克服图像噪声。

1 方法

本研究工作以CTA图像数据为实验材料，分三步对冠脉钙化进行量化：冠脉血管树重建、冠脉钙化分割、钙化分计算。

1.1 CTA数据获取

实验所用CTA数据来源于中国人民解放军254医院的GEⒸLightSpeedTM64排螺旋增强CT的心脏影像，其横向分辨率为512像素×512像素，纵向空间层间距为0.625 mm。

1.2 冠脉血管树重建

冠状动脉的分割与重建是冠脉钙化探测和量化的基础。在完成升主动脉分割［15］的基础上，本研究使用动态轮廓——Snake模型［16-17］来拟合升主动脉横截面轮廓，以实现对冠脉开口（轮廓上的异常突起）的识别，进而确定冠脉位置，使用基于自适应阈值［18］的区域生长算法完成对包含钙化在内的冠状动脉的分割，算法框架如图1所示。

1.3 冠状动脉钙化斑块分割

冠脉钙化分割总共包含3个步骤：血管区域划分、钙化斑块分割以及钙化分计算。

1.3.1 基于FCM聚类算法的血管区域划分

为克服造影剂引入的噪声，提高冠脉钙化量化的精度，本研究首先对冠脉血管内的像素点构建特征向量，然后根据血管形态计算自适应聚类数，进而采用FCM算法对冠脉血管进行区域划分。

1973年，Bezdek 等提出了 FCM 算法［19］，其代价函数与迭代公式分别为

图1 全自动冠脉识别分割算法框架。（a）原始的CTA图像；（b）基于Hough变换的升主动脉位置识别；（c）升主动脉分割结果；（d）冠脉开口识别结果（蓝色部分为冠脉开口）；（e）冠脉分割结果Fig.1 Algorithm frame for automatic coronary segmentation.（a）OriginalCTA image；（b）Location of aorta using Hough circle detection；（c）Segmented result of ascending aorta；（d）Recognition ofcoronary artery（volumes in blue）；（e）Segmenting resultofcoronary artery based on region growing

式中：uij表示第j个样本对第i个类的隶属度，m为聚类数，n为样本数，Zi为第i个的模糊类的聚类中心，dij表示第i个聚类中心与第j个样本之间的欧氏距离；α∈［1，∞），为算法的柔性参数，控制聚类的模糊程度。

为保证m值满足钙化探测需求，在研究工作中设置聚类数，其中r为CTA图像的像素间距，s为冠脉血管区域骨架化［20］后的像素数目。以X为聚类向量，取α＝2，采用FCM算法对冠脉血管内体素点进行聚类，选取样本最大隶属度对应的模糊组号作为该样本的最终隶属类，则冠脉区域划分结果如图2所示（CPR视图）。可以看出，冠脉区域内灰度分布相近的体素点被聚在同一区域，而噪声点则与相邻的体素点被聚在同一区域。

图2 基于FCM聚类算法的冠脉血管区域划分结果（CPR视图，蓝线为冠脉血管边界，红线为候选区域边界）Fig.2 Volume-divide result based on FCM algorithm in CPR view（blue lines indicate the boundary of the vessel and the red lines indicate the boundaries of volumes）

1.3.2 基于自适应阈值的钙化斑块分割

在本研究工作中，使用一维高斯函数对直方图峰值之后的部分进行拟合，得到最佳拟合参数μ与σ，其中一维高斯函数公式为

如图3所示，在本研究工作中，选取经验阈值T＝μ+1.5σ作为自适应阈值，对上述已经划分好的区域进行筛选，像素平均CT值大于阈值T的区域将会被认为是钙化斑块。图4分别显示了在三维视图和CPR视图下冠脉钙化探测的结果。

图3 自适应阈值选择方法Fig.3 Threshold selecting method

1.3.3 钙化分计算

在本研究工作中，使用钙化体积和Agatston钙化分来量化冠状动脉钙化，其中钙化体积计算方式如下：

图4 钙化探测结果。（a）三维视图下冠脉钙化探测结果（箭头指示了小型钙化斑块）；（b）CPR视图下冠脉钙化探测结果（红线指示了冠脉血管边界，蓝线指示了钙化斑块边界）Fig.4 Calcium detection result.（a）Calcium detection result in 3D view（the arrows indicate small calcium plaques extracted）；（b）Calcium detection result in CPR view（Red lines indicate the boundary ofthe vesseland the blue lines indicate the boundaries of calcium plaques extracted）

式中，a0是单个像素的实际体积，nv是钙化像素点的数目。

Agatston钙化分的计算方式［21］如下：

式中，a0是单个像素的实际体积，nv是钙化像素点的数目，f（v）为像素CT值对应的权重。

1.4 算法验证与统计分析流程

在实验阶段，采用50组CTA数据对本算法的钙化量化精度、阈值选择鲁棒性以及抗噪性能进行了测试。根据上述流程，在普通电脑（CPU IntelⒸ CoreTMi3-3220，3.30 GHz，RAM 8.00 GB）上采用Matlab 2014b进行编程实验。在实验中，采用的标准结果由3位经验丰富的放射科医师人工标记：对于冠脉区域内的一个像素点，当有两位或两位以上医师将其标记为钙化点时，该点将作为钙化点存入标准结果，否则为非钙化点。在50组实验数据中，30组被标记为含有钙化斑块，18组被标记为无钙化，其余2组由于存在严重的伪影和造影剂灰度过低的问题无法用于实验。

为了对本算法的灵敏度、特异性和精确性进行定量评价，将算法的分割结果与标准结果进行对比，计算了算法的真阳性值（true positives，TP）、真阴性值（true-negatives，TN）、假阳性值（false-positives，FP）与假阴性（false negatives，FN）等统计学参数［8］，并据此评价算法的灵敏度（sensitivity）Se、特异性（specificity）Sp与精确度（accuracy）Ac，其表达式分别为

在实验验证阶段，首先对本方法的钙化探测精确度和抗噪性能进行了统计学分析和实验验证（见本文第2.1节），然后对钙化分的计算精度进行了分析和实验比较（见本文第2.2节），最后与DBTD算法比较，对本方法阈值选择的鲁棒性进行了实验分析（见本文第2.3节）。

2 结果

2.1 算法精确度及抗噪性能实验结果

使用本算法对30组被放射科医师标记为含有钙化斑块的CTA实验数据进行灵敏度、特异性和精确度测试，三者对应的平均值分别达到了89.5%、98.6%与98.0%，具体结果如图5所示，而被标记为无钙化的18组数据也未出现被算法误标记钙化的情况。

在聚类过程中，由于具有较高CT值的噪声点与周围CT值较低的像素点距离较近而被划分到同一区域，于是区域的平均CT值会被这些像素点拉低，从而不会被算法标记为钙化。因此，本算法几乎不会误标记钙化斑块。图6所示为本算法在噪声数据中的测试结果，表明具有较好的抗噪性能，能够有效克服由造影剂引入的图像噪声。

图5 含钙化的30组CTA数据中算法的灵敏度、特异性与精确度量化结果Fig.5 Sensitivity，specificity and accuracy of the method in the calcified 30 data

图6 同阈值下本算法与阈值分割算法的抗噪性能比较。（a）放射科医师标记的标准钙化分割结果；（b）基于本算法的冠脉钙化分割结果；（c）基于阈值分割的冠脉钙化分割结果（图中箭头指示了由于造影剂使用引入的图像噪声）Fig.6 Anti-noise ability comparison of our method and direct threshold segmentation using the same threshold.（a）Standard resultmarked by the radiologist；（b）Calcium detection result using our method；（c）Direct threshold segmentation result（the arrows indicate the image noise voxels brought in by the contrast medium）

2.2 钙化分计算实验结果

钙化分是预测冠脉疾病的有效参数。在实验中，对比了本方法与 DBTD算法在钙化体积与Agatston钙化分量化上的精度。如图7所示，虽然使用DBTD算法得到的量化结果在一些数据集中与标准结果相近（如1、3、10等数据集），但是如果数据集的直方图下降趋势不明显或者具有明显的图像噪声，则该算法的鲁棒性将不可避免地受到影响。很明显，本算法在冠脉钙化中具有相对较好的精确度。进一步计算DBTD算法与本算法在钙化体积和Agatston钙化分量化上与标准结果的皮尔逊相关系数，如表1所示。显然，本算法的计算结果与标准结果具有更好的线性关系。

图7 本算法与DBTD算法的钙化量化效果比较。（a）钙化体积量化结果对比；（b）Agatston钙化分量化对比Fig.7 Calcium quantification compare of our method and DBTD method。（a）Quantitative comparison of calcium volume；（b）Quantitative comparison of Agatston score

表1 DBTD算法、本算法与标准结果的皮尔逊相关系数对比Tab.1 Pearson correlation coefficient quantitative comparison of standard result and method-mark result___

2.3 鲁棒性阈值选择实验结果

本研究采用高斯函数拟合冠脉区域直方图的方式，构建得到了更加稳健有效的分割阈值。在直方图波动强烈的情况下，对比本方法与文献［14］中提出的同类型的DBTD算法在分割阈值选择上的鲁棒性，如图8所示。

图8 DBTD算法与本算法在直方图波动明显时的冠脉钙化分割效果对比。（a）放射科医师标记的标准钙化提取结果；（b）基于DBTD算法的阈值选择；（c）基于DBTD算法的钙化分割结果；（d）基于本算法的阈值选择；（c）基于本算法的钙化分割结果Fig.8 Calcium extraction result comparison of our method and DBTD for CTA data with a complex histogram.（a）Standard calcium extraction result marked manually by the radiologist；（b）Threshold selection result using DBTD；（c）Corresponding calcium extraction result using DBTD；（d）Threshold selection result using our method；（e）Corresponding calcium extraction result using our method

从图8中可以看出，由于平滑后直方图曲线的下降趋势变缓，造成了DBTD算法错误的阈值选择，而本研究基于高斯函数拟合的阈值选择方法则不会受到影响，依然得到了比较合适的阈值。这是因为在DBTD算法中，算法通过计算冠脉直方图（最高峰值以后的部分）的一阶微分，并检验其是否超过经验阈值进行阈值选择；当直方图相对简单、平滑且下降趋势明显时，这种阈值选择方式能够得到合适的阈值，而当直方图相对复杂、下降趋势不明显或部分呈现停顿趋势时，DBTD算法极可能会造成错误的阈值选择。与此相比，本研究基于高斯函数拟合构造的自适应阈值兼顾了直方图变化的整体趋势，具有更好的鲁棒性。

3 讨论

本研究提出了一种基于FCM聚类与自适应阈值的冠脉钙化分割与量化方法，实现了全自动的冠脉钙化斑块分割与钙化分量化。该方法在保证较高灵敏度和特异性的同时，有效克服了由造影剂使用引入的图像噪声。该研究成果有助于促进计算机辅助手段在冠心病诊断领域的应用，为放射科医生提供可靠的CT数据信息，降低了由于生理疲劳等人为因素造成的误诊率。

从实验结果来看，该钙化量化方法的灵敏度和特异性分别达到了89.5%与98.6%，对应的钙化体积和Agatston钙化分与标准结果的皮尔逊系数也分别达到了0.974与0.975，而且钙化斑块分割结果有着较好的抗噪能力。该方法的高灵敏度与良好的抗噪能力主要基于以下两点实现：一是基于FCM算法的血管区域划分，二是基于高斯拟合的自适应分割阈值选取。

本研究在阈值分割之前引入FCM模糊聚类算法，将冠脉血管划分为CT值分布相似的区域。在聚类过程中，具有较高CT值的噪声点由于与周围CT值较低的像素点距离较近而被划分到同一区域，于是区域的平均CT值会被这些像素点拉低，从而不会将噪声点误标记为钙化，因此该算法具有较强的抗噪声能力。另外，本研究使用高斯函数拟合冠脉区域直方图，进而根据拟合参数构建得到了更加稳健有效的钙化分割阈值，克服了直方图波动对阈值选择鲁棒性造成的影响。

冠脉钙化量化是计算机辅助手段在冠心病诊断方面的重要应用。Eilot等提出的基于模型的钙化量化方法［12］与Wolterink等提出的基于模式识别的钙化识别方法［13］虽然能实现对冠脉钙化的分割和识别，但过程复杂，运算消耗大，不适合实际应用；Teßmann 等提出的 DBTD 算法［14］能够快速分割冠脉钙化，但其阈值选择的鲁棒性较差；本方法将FCM聚类算法应用于血管区域划分，并通过构造自适应阈值的方式，大大提高了阈值选择的鲁棒性以及算法的灵敏度和抗噪性能。

冠脉钙化的分割与量化不仅能为冠脉病变程度提供有效的评估依据，也是进一步检测冠脉高危病变的基础。作为相对稳定的冠脉病变，冠脉钙化的分割与量化并不是笔者研究工作的终点，今后会以此为基础，继续研究冠脉高危病变的检测和分类，最终构建完整的冠脉疾病检测和评价体系。

4 结论

本研究提出了一种基于聚类与自适应阈值的冠脉钙化分割与量化方法。

该算法能够实现自动钙化分计算，且其提供的钙化斑块的位置、大小、形状等信息对冠脉疾病的介入式治疗有着积极有效的意义。同时，对钙化斑块的分割也是进一步分类更危险病变（如冠脉粥样硬化等疾病）的基础。迄今为止，本算法在实验平台上处理一套CTA数据的时间为10～15 s，如果能够将其集成为专业软件并用于CT扫描仪的工作台，则时间会进一步缩短，因此有着潜在的商业应用价值。

［1］ Budoff MJ，Shaw LJ，Liu ST，et al.Long-term prognosis associated with coronary calcification：observations from a registry of 25，253 patients.［J］.Journal of the American College of Cardiology，2007，49（18）：1860-1870.

［2］ Polonsky TS，Mcclelland RL，Jorgensen NW，et al.Coronary artery calcium score and risk classification for coronary heart disease prediction. ［J］.Jama the Journal of the American Medical Association，2010，303（16）：1610-1616.

［3］ 乔爱科，侯映映，侯阳.冠状动脉狭窄几何构型对血流储备分数影响的有限元分析［J］.中国生物医学工程学报，2015，34（2）：198-203.

［4］ Sunkara N，Wong N D，Malik S.Role of coronary artery calcium in cardiovascular risk assessment. ［J］. Expert Review of Cardiovascular Therapy，2014，12（1）：87-94.

［5］ Seitun S，Maffei E，Martini C，et al.Calcium score and coronary plaque［M］//Clinical Applications of Cardiac CT.Milan：Springer Milan，2012：115-137.

［6］ Blaha MJ，Silverman MG，Budoff MJ.Is there a role for coronary artery calcium scoring for management of asymptomatic patients at risk for coronary artery disease？：Clinical risk scores are not sufficient to define primary prevention treatment strategies among asymptomatic patients.［J］.Circulation Cardiovascular Imaging，2014，7（2）：398-408.

［7］ 何玉娜，蔺嫦燕.冠脉支架内再狭窄的血流动力学研究进展［J］. 中国生物医学工程学报，2015，34（3）：354-359.

［8］ Gouw DJJMD. Automatic coronary calcium scoring using computed tomography［J］.2015.

［9］ Glodny B，Helmel B，Trieb T，et al.A method for calcium quantification by means of CT coronary angiography using 64-multidetector CT：Very high correlation with Agatston and volume scores.［J］.European Radiology，2009，19（7）：1661-1668.

［10］ Schuhbaeck A，Otaki Y，Achenbach S，et al.Coronary calcium scoring from contrast coronary CT angiography using a semiautomated standardized method ［J］. Journal of Cardiovascular Computed Tomography，2013，61（10）：E914.

［11］ Otton JM，Lønborg JT，Boshell D，et al.A method for coronary artery calcium scoring using contrast-enhanced computed tomography ［J］. Journal of Cardiovascular Computed Tomography，2012，6（1）：37-44.

［12］ Eilot D，Goldenberg R.Fully automatic model-based calcium segmentation and scoring in coronary CT angiography.［J］.International Journal of Computer Assisted Radiology＆ Surgery，2014，9（4）：595-608.

［13］ Wolterink JM，Leiner T，de Vos BD，et al.Automatic coronary artery calcium scoring in cardiac CT angiography using paired convolutional neural networks.［J］.Medical Image Analysis，2015，9349：589-596.

［14］ Teβmann M，Vega-Higuera F，Bischoff B，et al.Automatic detection and quantification of coronary calcium on 3D CT angiography data ［J］. Computer Science-Research and Development，2011，26（26）：117-124.

［15］ 杨栋.面向CTA图像的冠脉血管分割算法研究和血管狭窄度分析［D］.杭州：浙江大学，2015.

［16］ Kass M，Witkin A，Terzopoulos D.Snakes：Active contour models［J］.International Journal of Computer Vision，1988，1（4）：321-331.

［17］ Cohen LD.On active contour models and balloons［J］.Cvgip Image Understanding，1991，53（2）：211-218.

［18］ 张朝霞.基于MSCT的感兴趣血管段CAG最佳视角和冠脉运动估计的研究［D］.天津：天津大学，2011.

［19］ James C.Bezdek†.Cluster validity with fuzzy sets［J］.Journal of Cybernetics，1973，3（3）：58-73.

［20］ Lee TC，Kashyap RL，Chu CN.Building skeleton models via 3-D medial surface axis thinning algorithms［J］.Cvgip Graphical Models＆ Image Processing，1994，56（6）：462-478.

［21］ Agatston AS，Janowitz WR，Hildner FJ，et al.Quantification of coronary artery calcium using ultrafast computed tomography［J］.Journal of the American College of Cardiology，1990，15（4）：827-832.

Segmentation and Scoring of Coronary Calcium in 3D CTA Data Based on FCM Algorithm and Self-Adapting Threshold Determination

Zhao Cong1Chen Xiaodong1∗Zhang Jiachen1Wang Yi1Jia Zhongwei2Chen Xiangzhi3

1（Key Laboratory of Opto-Electronic Information Science and Technology of Ministry of Education，School of Precision Instrument and Opto-Electronics Engineering，Tianjin University，Tianjin300072，China）

2（Department of Cardiovascular Medicine，People′s Liberation Army254Hospital，Tianjin300142，China）

3（Department of Radiology，People′s Liberation Army254Hospital，Tianjin300142，China）

Aiming at the problems of image noise and impressionable threshold in coronary calcium quantification with 3D CTA data，a new method based on fuzzy C-means（FCM）clustering algorithm and selfadapting threshold determination was proposed for automatically segmenting and quantifying the coronary calcium.Firstly，feature vectors are constructed for every voxel in the coronary artery，which contains spatial coordinates and CT value information of the voxel，and a clustering algorithm combined with a self-adapted group number referring to vessel skeleton number is used to divide the coronary artery into different candidate volumes；secondly，a robust threshold determination algorithm based on the histogram is used to extract calcium plaques among those candidate volumes acquired；finally，calcium volume and Agatston score are calculated.Result shows that the method proposed in this study has relatively high sensitivity of 89.5%and specificity of 98.6%in calcium detection in 30 coronary CTA data.The calcium volume and Agatston score calculated automatically show a high correlation with the standard result.The corresponding Pearson correlation coefficient being up to 0.974 and 0.975，respectively，much higher than the 0.523 and 0.501 that calculated by the derivative-based threshold determination（DBTD）method.Experimental results show that this method can be used for coronary calcium segmentation and quantification，and has the characteristics of full automation，robustness and noise immunity.

CT image processing；coronary calcium；segmentation and quantification；fuzzy C-means（FCM）clustering algorithm；calcium score calculation

R318 文献标志码：A 文章编号：0258-8021（2017）05-0550-07

10.3969 /j.issn.0258-8021.2017.05.006

2016-10-13，，录用日期：2017-06-06

∗通信作者（Corresponding author），E-mail：xdchen@tju.edu.cn

（致谢：感谢中国人民解放军第254医院浦奎、王星医生对本论文研究工作提供的帮助）