基于多特征支持向量机和弹性区域生长的膝软骨自动分割

引用本文

王品, 何璇, 吕洋, 李勇明, 邱明国, 刘书君. 基于多特征支持向量机和弹性区域生长的膝软骨自动分割. 2016, 46(5): 1688-1696
WANG Pin, HE Xuan, LYU Yang, LI Yong-ming, QIU Ming-guo, LIU Shu-jun. Automatic segmentation of articular cartilages using multi-feature SVM and elastic region growing. Journal of Jilin University Engineering and Technology Edition, 2016, 46(5): 1688-1696 复制到剪切板

Permissions

基于多特征支持向量机和弹性区域生长的膝软骨自动分割

王品¹, 何璇¹, 吕洋¹, 李勇明^1,², 邱明国², 刘书君¹

1.重庆大学通信工程学院,重庆 400044

2.第三军医大学生物医学工程学院,重庆 400038

作者简介:王品(1979-),女,副教授,博士.研究方向:医学图像处理.E-mail:wangpin@cqu.edu.cn

基金:国家自然科学基金项目(61108086,11304382); 中国博士后科学基金项目(2013M532153); 重庆市自然科学基金项目(CSTC 2012jjA40015); 重庆市科技攻关计划项目(cstc2012gg-yyjs0572); 中央高校基本科研业务费项目(CDJZR12160011, CDJZR13160008, CDJZR155507).

摘要

为了从膝关节磁共振图像(MRI)中分割出膝软骨,提出一种基于多特征支持向量机(SVM)边缘定位和弹性区域生长的自动分割算法。首先,采用自适应Canny边缘检测算法提取图像主要边缘;再对边缘提取多个图像特征,结合SVM算法对边缘进行分类,完成软骨边缘定位;然后,在软骨边缘的基础上进行种子点及软骨像素区域的选择;之后基于选择的结果采用弹性区域生长进行初步软骨分割;最后,基于先验知识和形态学获得最终膝软骨分割结果。实验结果表明:该算法能够准确、快速地自动分割出膝关节MRI中不同的膝软骨,其中股软骨、胫软骨、髌软骨的平均评价重要指标(DSC)分别可达0.8543、0.8280、0.8703,与手工分割结果具有较高的一致性。

关键词: 信息处理技术; 多特征支持向量机; 弹性区域生长; 膝软骨自动分割

中图分类号:TN911.73;R445.2 文献标志码:A 文章编号:1671-5497(2016)05-1688-09

Automatic segmentation of articular cartilages using multi-feature SVM and elastic region growing

WANG Pin¹, HE Xuan¹, LYU Yang¹, LI Yong-ming^1,², QIU Ming-guo², LIU Shu-jun¹

1.College of Communication Engineering,Chongqing University, Chongqing 400044,China

2.College of Biomedical Engineering, The Third Military Medical University, Chongqing 400038,China

Abstract

We present an algorithm that automatically and accurately segments the cartilages from Magnetic Resonance Images (MRI) of knees using multi-feature Support Vector Machine (SVM) and elastic region growing. First, adaptive canny edge detection is used to extract the main edges of the images. Second, the edge features are extracted and the edges are classified with SVM to complete the location of the cartilage edges. Then, based on the cartilage edges, the seeds and regions of the cartilage are chosen, and the preliminary segmentation of knee cartilages is achieved using the elastic region growing. Finally, the priori knowledge and morphology are used to improve the preliminary segmentation results. Experimental results show that the proposed method can accurately and rapidly segment the knee cartilages with high consistency with manual segmentation results. The proposed method obtains average DSC values of 0.8543, 0.8280 and 0.8703 for the femoral, tibial and patellar cartilages.

Keyword: information processing technology; multi-feature support vector machine; elastic region growing; automatic segmentation of knee cartilages

Show Figures

0 引言

骨关节炎(OA)是最常见的膝关节失调疾病, 危害很大。近年来发现, OA的患者已经出现低龄化趋势, 目前年龄最小的只有25岁, 因此需要监测其退化程度并实施早期诊断, 从而使其得以有效治疗。研究表明^[1], OA的发作和退化等级与膝关节软骨(AC)密切相关。由于AC的损害程度反映疾病的严重程度, 因此临床常用AC的损害程度作为诊断和评价OA退化程度的重要指标^[2]。磁共振图像(MRI)具有无创性、无辐射性以及良好组织分辨率等性质, 已经被证明是监测AC损害程度(即形态及生理变化)的最适宜的成像检测手段^{[3, 4]}。

基于膝关节MRI进行AC的准确分割, 对AC的损害程度的定量测量至关重要。只有当软骨能准确分割出来, 其厚度、体积、表面积等参数才能被准确计算, 从而辅助医生进行膝关节疾病诊断, 为制定手术方案提供有力的数据基础。目前, 临床上通常采用全手动或半自动的AC分割方法^{[5, 6, 7, 8]}, 这些方法费时费力、效率较低, 因此设计一种面向AC精准、快速的全自动分割方法显得尤为迫切。然而, 由于软骨的形态细长且薄, 软骨与其他组织之间对比度较低, MRI中存在伪影及灰度不均等原因, 因此从MRI中自动分割出软骨是一个重要难题。近几年, 许多学者对全自动AC分割算法进行了探索研究。Lee等^[9]使用图谱和局部结构分析方法提取关节软骨; González等^[10]提出了一种基于主动轮廓模型(ASM)与局部二进制模式(LBP)相结合的半自动分割方法; Prasoon等^{[11, 12]}提出了3个二维卷积神经网络(CNNs)组成的新型像素分类整合系统, 这3个网络分别对应三维图像的 $xy 、 yz$ 和 $\begin{matrix} zx \end{matrix}$ 平面。目前, 国内对于AC自动分割的研究较少。庞剑飞等^{[13, 14]}使用边缘定位与支持向量机(SVM)相结合的方法自动分割软骨。

本文分析了膝关节MRI中AC的特点, 提出了一种考虑空间相关性关系的多特征提取及模式分类法来自动识别骨-软骨边缘, 并在边缘定位的基础上, 采用弹性自动区域生长算法实现软骨的自由分割。该算法充分利用了软骨解剖学信息, 通过边缘分类来实现精确定位, 并结合弹性区域生长获得精细的分割效果和可接受的时间代价。

1 图像数据来源

本文采用西南医院提供的5名无关节病史的健康成年男性(平均年龄23.3岁, 年龄在20~25岁)右膝关节矢状位MRI图像作为研究对象, 图像使用1.5T西门子扫描仪得到, 采用T2加权脂肪抑制序列(回波时间为4.42 ms, 重复时间为1363 ms, 偏转角为60° , 切片厚度为2.5 mm, 视场为160 mm, 矩阵尺寸为384× 384, 层数为25左右)。T2加权脂肪抑制序列是临床上进行关节软骨检查的常用序列, 能较为清晰地显示髌软骨、股软骨和胫软骨。图1为某受检者的膝关节图像。

	Figure Option View Download New Window
	图1 膝关节MRI原始图像Fig.1 MRI of knee

从图1可以看出, 股软骨、胫软骨和髌软骨的一侧分别与股骨、胫骨和髌骨相连, 两者之间有较为明显的灰度差异, 但另一侧股软骨与髌软骨相连, 股软骨和胫软骨与半月板、肌肉及脂肪等组织相连, 其灰度相近, 且软骨的形态细长且薄, 在分割时很容易受到周围复杂纹理的影响。观察MRI图像可以发现, 股骨、胫骨和髌骨的灰度与周围灰度差异较大, 因此, 首先提取出骨-软骨边缘, 然后, 在骨-软骨边缘基础上分割出软骨。然而, 由于骨骼内部有噪声干扰, 骨边缘灰度不均, 股软骨与髌软骨、胫软骨分别相连, 此外, 肌肉中的某些灰度也明显高于周围组织, 综合种种原因, 导致单纯地采用阈值或Canny检测的方法很难直接得到骨-软骨边缘。

2 多特征支持向量机和弹性区域生长的膝关节软骨自动分割算法

2.1 算法概述

本文算法的总体流程图如图2所示。由于软骨及骨骼周围各组织形态及纹理的影响, 使得软骨与非软骨之间存在边界模糊的情况。本文提出在骨-软骨边缘定位的基础上, 采用改进的弹性区域生长算法对膝关节软骨进行精确分割, 该方法主要分为两个步骤:①基于多特征实现膝关节MRI图像中各骨-软骨边缘定位; ②基于提取的骨-软骨边缘定位, 通过自动选取种子点的弹性区域生长算法对软骨进行精确分割。

	Figure Option View Download New Window
	图2 膝软骨自动分割算法总体流程图Fig.2 Flow chart of the automatic segmentation of knee cartilages

2.2 骨-软骨边缘提取

骨-软骨边缘的提取主要分为:自适应Canny边缘检测、边缘特征提取、边缘识别。

2.2.1 自适应Canny边缘检测Canny算子使用两个阈值来分别检测强边缘和弱边缘, 而且仅当弱边缘与强边缘相连时, 弱边缘才会包含在输出中。因此, 该方法不容易受噪声干扰, 能够检测到弱边缘。但是, 传统Canny算子的高、低阈值 $\begin{matrix} T_{h} \end{matrix}$ 和 $\begin{matrix} T_{l} \end{matrix}$ 需手动设置, 而膝关节MRI图像的切片层数在25层左右, 手动设置 $\begin{matrix} T_{h} \end{matrix}$ 和 $\begin{matrix} T_{l} \end{matrix}$ 不仅耗时, 而且单纯地以人的感官选取会存在主观性强、阈值普适性差等问题, 因此本研究采用自适应Canny边缘检测算法^[15], 动态选取合适的 $\begin{matrix} T_{h} \end{matrix}$ 和 $\begin{matrix} T_{l} \end{matrix}$ 以提高后续的边缘检测效率和普适性。本文采用如下改进的自适应Canny边缘检测算法:

Step1 对图像进行高斯平滑, 计算平滑后图像的梯度幅值图像, 并对梯度幅值图像应用非最大抑制得到新的梯度幅值图像 $\begin{matrix} Imgra d_{NMS} 。 \end{matrix}$

Step2 对 $\begin{matrix} Imgra d_{NMS} \end{matrix}$ 采用迭代算法^[16]得到全局阈值 $\begin{matrix} T_{grad}, \end{matrix}$ 并计算 $\begin{matrix} T_{grad} \end{matrix}$ 左右两个区域的均值 $\begin{matrix} μ_{l} \end{matrix}$ 和 $\begin{matrix} μ_{h}, \end{matrix}$ 以及方差 $\begin{matrix} σ_{l} \end{matrix}$ 和 $\begin{matrix} σ_{h}, \end{matrix}$ 动态地设定高、低阈值为:

$\begin{matrix} \{\begin{matrix} T_{h} = μ_{h} + λ_{h} \times σ_{h} \\ T_{l} = μ_{l} + λ_{l} \times σ_{l} \end{matrix} (1) \end{matrix}$

最后, 根据边缘的数目动态地确定系数 $\begin{matrix} λ_{h} \end{matrix}$ 和 $\begin{matrix} λ_{l} \end{matrix}$ 的值, 边缘过多则提高 $\begin{matrix} λ_{h} \end{matrix}$ 和 $\begin{matrix} λ_{l}, \end{matrix}$ 边缘过少则降低 $\begin{matrix} λ_{h} \end{matrix}$ 和 $\begin{matrix} λ_{l}, \end{matrix}$ 最终得到Canny算子的 $\begin{matrix} T_{h} \end{matrix}$ 和 $\begin{matrix} T_{l} 。 \end{matrix}$

图3为某受试者MRI图像的第3、7、15层切片。采用上述改进的自适应Canny边缘检测算法动态设定高、低阈值对图像进行边缘检测, 初始结果如图3(d)(e)(f)所示。根据解剖学知识, 结合髌骨位置及形态的固定性等特点, 对髌骨边缘的角点进行断开, 断开边缘中的分叉点如图3(g)(h)(i)所示。

	Figure Option View Download New Window
	图3 自适应Canny边缘检测结果Fig.3 Results of adaptive Canny edge detection

2.2.2 边缘特征提取

由图3可以看出, 采用上述自适应Canny边缘检测算法, 除了检测出骨边缘外还会检测出肌肉、组织等非骨边缘, 因此本文采用支持向量机(SVM)对检测出的边缘进行真伪边缘的分类, 从而识别出股骨边缘、胫骨边缘和髌骨边缘。使用SVM对边缘进行分类时首先需要进行边缘的特征提取, 而提取的特征直接影响分类的精度。根据MRI图像的特点及解剖学知识, 本文提取出边缘的局部信息、几何学信息以及邻近体素之间的信息, 共有20个特征, 如表1所示。分析表1中各特征的意义可知:特征1、4及6~11代表了图像的局部信息, 而特征2、3、5代表了邻近体素之间的信息, 特征12~20代表了图像的几何学信息, 因此提取的特征较为丰富, 比较全面地描述了边缘的特性。

表1 边缘特征 Table 1 Features of edge

编号	特征属性	特征公式	特征描述
1	线上点的平均灰度值	$\begin{matrix} \bar{f} = [\sum_{x, y \subset L} f (x, y)] / n \end{matrix}$	$\begin{matrix} L \end{matrix}$ 为边缘线; $\begin{matrix} f \end{matrix}$ 为原图的灰度值; $\begin{matrix} n \end{matrix}$ 是边缘线上的像素点个数。
2	线上每个点5邻域的平均灰度值	$\begin{matrix} \bar{f} = [\sum_{x, y \subset L_{5}} f (x, y)] / n_{5} \end{matrix}$	$\begin{matrix} L_{5} \end{matrix}$ 为边缘线上像素点周围5邻域范围; $\begin{matrix} n_{5} \end{matrix}$ 为一条边缘线上所有像素点5邻域范围内像素点总数。
3	线上每个点5邻域方差的均值	$\begin{matrix} \bar{v} = [\sum_{x, y \subset L} va r_{5} (f (x, y))] / n \end{matrix}$	var₅为线上像素点5邻域范围内灰度值的方差。
4	线上灰度值的方差	$\begin{matrix} v = var (f (x, y)) x, y \subset L \end{matrix}$	表示边缘线上像素点之间的灰度值方差。
5	线上5邻域灰度值的方差	$\begin{matrix} v_{5} = var (f (x, y)) x, y \subset L_{5} \end{matrix}$	表示边缘线上像素点5邻域内像素间的灰度值方差。
6	每条边缘线的长度	$\begin{matrix} n \end{matrix}$	$\begin{matrix} n \end{matrix}$ 为每条边缘线上像素点总和。
7	线上点的梯度幅值的均值	$\begin{matrix} \bar{g} = (\sum_{x, y \subset L} grad (f (x, y))) / n \end{matrix}$	$\begin{matrix} grad \end{matrix}$ 为边缘线上像素点的梯度幅值。
8	平均的 $\begin{matrix} X \end{matrix}$ 坐标	$\begin{matrix} \bar{x} = (\sum_{x, y \subset L} x) / n \end{matrix}$	$\begin{matrix} x \end{matrix}$ 为边缘线上像素点的 $\begin{matrix} X \end{matrix}$ 坐标。
9	平均的 $\begin{matrix} Y \end{matrix}$ 坐标	$\begin{matrix} \bar{y} = (\sum_{x, y \subset L} y) / n \end{matrix}$	$\begin{matrix} y \end{matrix}$ 为边缘线上像素点的 $\begin{matrix} Y \end{matrix}$ 坐标。
10	$\begin{matrix} X \end{matrix}$ 坐标的方差	$\begin{matrix} v_{x} = var (x) x, y \subset L \end{matrix}$	表示边缘线上像素间 $\begin{matrix} X \end{matrix}$ 坐标方差。
11	$\begin{matrix} Y \end{matrix}$ 坐标的方差	$\begin{matrix} v_{y} = var (y) x, y \subset L \end{matrix}$	表示边缘线上像素间 $\begin{matrix} Y \end{matrix}$ 坐标方差。
12	水平投影	$\begin{matrix} w = \max (x) - \min (x) \end{matrix}$	表示边缘线在 $\begin{matrix} x \end{matrix}$ 轴水平投影的大小。
13	垂直投影	$\begin{matrix} h = \max (y) - \min (x) \end{matrix}$	表示边缘线在 $\begin{matrix} y \end{matrix}$ 轴垂直投影的大小。
14	左侧边缘像素数	$\begin{matrix} numleft = \sum_{x, y \subset L} (A (y, :) < A (y, x)) \end{matrix}$	表示在边缘线左侧的边缘点数目, $\begin{matrix} A \end{matrix}$ 为所有边缘的图像矩阵。
15	右侧边缘像素数	$\begin{matrix} numright = \sum_{x, y \subset L} (A (y, :) > A (y, x)) \end{matrix}$	表示在边缘线右侧的边缘点数目。
16	开口方向	$\begin{matrix} opendirect = \max (T_{up}, T_{down}, T_{colliner}) \end{matrix}$	T_up、T_down、T_colliner分别为开口向上、向下、共线的比率。
17	旋转角度	$\begin{matrix} angle = abs (angl e_{1} - angl e_{end}) \end{matrix}$	$\begin{matrix} angl e_{1} 、 angl e_{end} \end{matrix}$ 分别为边缘的起始点和终点与水平方向夹角。
18	旋转方向	$\begin{matrix} angdirect = \{\begin{matrix} 1, & angle < 180 ° \\ - 1, & angle > 180 ° \end{matrix} \end{matrix}$	表示边缘线的旋转方向, 与水平方向夹角大于180° 为-1; 小于180° 为1。
19	水平像素数	$\begin{matrix} numhori = Freeman (1) + Freeman (5) \end{matrix}$	$\begin{matrix} Freeman (1) 、 Freeman (5) \end{matrix}$ 分别为Freeman编码中方向为0和4的点的数目。
20	左侧边缘数	$\begin{matrix} ledgenum = \sum ({\bar{x}}_{A} < \bar{x}) \end{matrix}$	表示边缘线左边存在的边缘数目, ${\bar{x}}_{A}$ 为边缘的平均 $\begin{matrix} x \end{matrix}$ 坐标。

表1 边缘特征 Table 1 Features of edge

2.2.3 边缘识别

SVM^[17]是一种监督机器学习方法, 将输入空间线性不可分的样本通过核函数映射到一个高维特征空间, 转换为在特征空间中求解线性约束下的凸二次规划问题, 目标是找到一个最优的分离超平面, 将所有的训练向量正确分类。假设一个训练集为 $\begin{matrix} \{(x_{i}, y_{i}) |i = 1, 2, \dots, N\}, \end{matrix}$ 其中, $\begin{matrix} x_{i} \end{matrix}$ 为特征向量, $\begin{matrix} y_{i} \end{matrix}$ 为 $\begin{matrix} x_{i} \end{matrix}$ 相应类的标识, $\begin{matrix} N \end{matrix}$ 为训练样本的数量, 则可以通过最大化式(2)计算出最优拉格朗日乘子 $\begin{matrix} λ_{i}, \end{matrix}$ 通过式(3)可以得到最优的超平面。本文采用的SVM分类模型如下:

$\begin{matrix} \begin{matrix} \{\begin{matrix} \max_{λ} (\sum_{i} λ_{i} - \frac{1}{2} \sum_{i, j} λ_{i} λ_{j} y_{i} y_{j} K (x_{i}, x_{j})) \\ 0 \leq λ_{i} \leq C, i = 1, 2, \dots, N \\ \overset{N}{\sum_{i = 1}} λ_{i} y_{i} = 0 \end{matrix} (2) \\ f (x) = \overset{N_{s}}{\sum_{i = 1}} λ_{i} y_{i} K (x_{i}, x) + b (3) \end{matrix} \end{matrix}$

式中: $\begin{matrix} x \end{matrix}$ 为测试数据; 偏差 $\begin{matrix} b \end{matrix}$ 是通过所有条件的均值计算出的; $\begin{matrix} x_{i} \end{matrix}$ 为具有最优解的向量参数; $N_{s}$ 为支持向量的数量; $λ$ 为拉格朗日系数; $\begin{matrix} C \end{matrix}$ 为惩罚因子。

核函数采用高斯径向基核函数(RBF):

$\begin{matrix} K (x_{i}, x_{j}) = \exp (- \frac{‖ x_{i} - x_{j} ‖^{2}}{σ^{2}}) (4) \end{matrix}$

本文的目标是分割出股骨边缘、胫骨边缘和髌骨边缘及非骨边缘, 属于多分类问题。本研究采用一对一方法^[18]来解决该多分类问题。一对一方法是在 $\begin{matrix} k \end{matrix}$ 类问题中进行两两组合, 构建 $\begin{matrix} k (k - 1) / 2 \end{matrix}$ 个不同的决策函数, 采用“ 投票机制” , 得出最终测试样本 $\begin{matrix} x \end{matrix}$ 所属的类 $\begin{matrix} i, \end{matrix}$ 当且仅当

$\begin{matrix} C_{x} = argmax \overset{k}{\sum_{j \neq i, j = 1}} sgn f_{ij} (x) (5) \end{matrix}$

式中: $\begin{matrix} j \end{matrix}$ 为分类器的编号; $\begin{matrix} C_{x} \end{matrix}$ 为测试数据 $\begin{matrix} x \end{matrix}$ 属于 $\begin{matrix} i \end{matrix}$ 类时的投票数; $\begin{matrix} sgn f_{ij} (x) \end{matrix}$ 为符号函数, 当 $\begin{matrix} f_{ij} (x) \end{matrix}$ 值为正数时, 函数值为1, 否则为0。

若多个类别的投票数相同, 则对投票相同的类别采用一对一的方法再次进行投票, 直到无投票数相同的情况为止。

本文需要将检测到的边缘分为股骨边缘、胫骨边缘、髌骨边缘及非骨边缘四类, 因此需要构造6个不同的两两分类器。将提取的特征分为三大类:代表图像局部信息的特征1、4及6~11, 代表邻近体素之间信息的特征2、3、5, 代表图像几何学信息的特征12~20。为了验证本文提取的特征的有效性, 分别选取不同的特征对5个人的128幅图像进行边缘识别, 特征的选取组合原则及识别准确率如表2所示。采用FV2虽然也能得到较高的识别率, 但某些图像中会出现虚假边缘。而采用改进后的特征FV3, 识别率达100%。

表2 采用不同特征的边缘识别率 Table 2 Edge recognition rate with different feature

为了保证识别的准确性, 识别完成后还需进行虚假边缘检测及后处理。

2.3 虚假边缘检测及边缘后处理

由于初始Canny边缘检测的骨边缘中只有部分像素点为骨-软骨边缘像素点, 因此需要对识别出的骨边缘进行后处理。根据先验知识, 考虑各骨骼之间的位置及距离关系去掉虚假边缘点。后处理步骤如下:

Step1 对于股骨和髌骨粘连的情况(见图3(b)), 初始Canny边缘检测会出现如图3第2列所示的情况, 为后续软骨分割的精确, 需要去掉非髌骨-软骨边缘的点。考虑位置的先验知识, 髌骨-软骨边缘应位于股骨边缘左侧, 设 $\begin{matrix} x_{p} 、 y_{p} \end{matrix}$ 和 $\begin{matrix} x_{f} 、 y_{f} \end{matrix}$ 分别代表髌骨和股骨边缘像素点的 $\begin{matrix} X \end{matrix}$ 和 $\begin{matrix} Y \end{matrix}$ 的坐标。因此, 需要去掉满足 $\begin{matrix} x_{p} = x_{f} \end{matrix}$ 且 $\begin{matrix} y_{p} > y_{f} \end{matrix}$ 条件的髌骨边缘像素点, 得到准确的髌骨-软骨边缘 $\begin{matrix} edg e_{p} 。 \end{matrix}$

Step2 根据解剖学知识, 股骨-软骨边缘应该在股骨头处, 根据几何学的知识, 考虑股骨重心与股骨边缘点及边缘点重心的位置关系, 去掉股骨右侧的非股骨-软骨接触面的像素点。设 $\begin{matrix} centr e_{f} \end{matrix}$ 为股骨的重心坐标, $\begin{matrix} centr e_{fl} \end{matrix}$ 为检测出的股骨-软骨边缘线上每隔 $\begin{matrix} p \end{matrix}$ 个点选取 $\begin{matrix} l \end{matrix}$ 长度的线段的重心坐标, 当两重心相连的线段与股骨边缘线相交时, 则断开交点, 并去掉交点后的点, 得到准确的股骨-软骨边缘 $\begin{matrix} edg e_{f} 。 \end{matrix}$

Step3 对于检测出的胫骨边缘, 采用Harris角点检测算法检测胫骨边缘中的角点。计算每个角点左右邻域及上下邻域灰度差 $\begin{matrix} f_{rl} \end{matrix}$ 和 $\begin{matrix} f_{ud} \end{matrix}$ 的灰度差, 若 $\begin{matrix} f_{rl} > T_{f} \end{matrix}$ 且 $\begin{matrix} f_{ud} < T_{f}, \end{matrix}$ 其中 $\begin{matrix} T_{f} \end{matrix}$ 为阈值, 则断开该角点, 并保留水平像素数较多的胫骨边缘线, 得到准确的胫骨-软骨边缘 $\begin{matrix} edg e_{t} 。 \end{matrix}$

经过SVM分类的骨边缘及后处理得到的骨-软骨边缘结果如图4所示, 其中1代表股骨-软骨边缘, 2代表胫骨-软骨边缘, 3代表髌骨-软骨边缘。

	Figure Option View Download New Window
	图4 提取出的骨-软骨接触面Fig.4 Detected bone-cartilage interface

2.4 结合骨-软骨接触面的区域生长

区域生长法^[19]的基本思想是:选取一组种子点, 将与种子预先定义的性质相似的那些邻域像素添加到每个种子上来形成这些生长区域。其主要优点是思想简单, 算法容易实现, 对于较均匀的具有相同特征的连通区域有较好的分割效果。区域生长算法的关键在于:①种子点的选取。传统区域生长法的种子点多为手动选择, 而膝关节MRI图像有25层左右, 需手动选择25次, 过于繁琐, 因此本文结合骨-软骨边缘, 自动选取种子点, 避免人为干扰, 保证了种子点选取的稳定性。②相似性准则的确定。传统的区域生长准则主要是基于区域灰度差, 将相似性准则定义为待生长的像素灰度值 $\begin{matrix} f (x, y) \end{matrix}$ 与整个图像区域的灰度均值 $\begin{matrix} m \end{matrix}$ 之差, 而软骨区域较小, 骨骼、背景等面积较大的区域又与软骨灰度差异较大, 使得全局灰度均值并不能表征软骨的灰度均值。因此, 本文采用生长变化的局部灰度均值 $\begin{matrix} m' \end{matrix}$ 来代替全局灰度均值 $\begin{matrix} m, \end{matrix}$ 动态地设定阈值 $\begin{matrix} T 。 \end{matrix}$ ③终止条件的确定。当不再有新的像素点进入生长区域时, 算法结束。本文采用的弹性区域生长的步骤如下:

Step1 结合骨-软骨边缘自动选取初始种子点。软骨厚度一般为1~6 mm, 相当于数字图像中的3~18个像素, 因此将候选初始种子点区域限定在骨-软骨周围20个像素内。分别计算股骨-软骨右侧, 胫骨软骨右侧及髌骨-软骨右侧20像素内的梯度幅值, 选取相应区域内梯度幅值及方差较小的最大连通区域 $\begin{matrix} fiel d_{i} (i = 1, 2, 3) \end{matrix}$ 作为初始种子点的候选区域, 在候选区域内随机选取 $\begin{matrix} m_{i} \end{matrix}$ 个种子点。虽然, 软骨区域灰度值较大且变化均匀, 但是在骨-软骨边缘周围依旧存在灰度变化较大区域, 为了确保初始种子点选取的全面性, 在骨-软骨边缘与 $\begin{matrix} fiel d_{i} \end{matrix}$ 之间的区域内随机选取 $\begin{matrix} n_{i} \end{matrix}$ 个初始种子点, 综合 $\begin{matrix} m_{i} \end{matrix}$ 和 $\begin{matrix} n_{i} \end{matrix}$ 得到最终的种子点 $\begin{matrix} see d_{i} (i = 1, 2, 3), \end{matrix}$ 其中, $\begin{matrix} i \end{matrix}$ 代表种子点的种类, 1代表股软骨种子点, 2代表胫软骨种子点, 3代表髌软骨种子点。

Step2 判断是否符合相似性准则。传统区域生长主要是将种子点邻域像素灰度值 $\begin{matrix} f (x, y) \end{matrix}$ 与整个图像区域的灰度均值 $\begin{matrix} m \end{matrix}$ 进行比较, 若两者之差小于 $\begin{matrix} T \end{matrix}$ 时, 则将该邻域像素点归并到一个区域。本文使用变化的局部区域均值 $\begin{matrix} m' \end{matrix}$ 代替全局均值 $\begin{matrix} m \end{matrix}$ , 定义相似性准则为:

$\begin{matrix} \{\begin{matrix} |f (x, y) - m'| < T \\ m' = \frac{1}{n'} \sum_{(x, y) \in R'} f (x, y) \end{matrix} (6) \end{matrix}$

式中: $\begin{matrix} R' \end{matrix}$ 为伴随区域生长变化的新局部区域; $\begin{matrix} n' \end{matrix}$ 为生长区域变化后的像素点总数; 阈值 $\begin{matrix} T \end{matrix}$ 的固定设置会影响分割的普适性, 因此本文将阈值 $\begin{matrix} T \end{matrix}$ 动态地设定为区域生长的标准差 $\begin{matrix} σ, \end{matrix}$ 表达式为:

$\begin{matrix} σ = \sqrt{\frac{1}{n'} \sum_{(x, y \in R')} {|f (x, y) - m'|}^{2}} (7) \end{matrix}$

Step3 判定Step2中满足条件的种子点是否在软骨候选区域内。由于软骨在骨-软骨边缘3~18个像素内, 为了减少算法冗余度, 设定区域生长的软骨候选区域分别为股骨-软骨边缘下方, 胫骨-软骨边缘上方及髌骨-软骨边缘右方20个像素内。计算Step2中满足条件的种子点与相应骨-软骨边缘的欧氏距离, 满足式(8)则保留种子点, 否则删除种子点。

$\begin{matrix} \min (\sqrt{{(x_{ij} - x'_{ij})}^{2} + {(y_{ij} - y'_{ij})}^{2}}) < T_{distance} (8) \end{matrix}$

式中: $\begin{matrix} (x_{ij}, y_{ij}) \end{matrix}$ 为相应软骨的第 $\begin{matrix} j \end{matrix}$ 个种子点的坐标; $\begin{matrix} (x'_{ij}, y'_{ij}) \end{matrix}$ 为相应骨-软骨边缘上的点的坐标; $\begin{matrix} T_{distance} \end{matrix}$ 为阈值。

Step4 判断是否满足终止条件。将满足式(8)的种子点作为新的种子点, 重复Step2和Step3, 继续寻找满足相似性准则的像素, 若没有新的点加入到区域生长中, 则终止。

2.5 软骨分割结果后处理

在本算法中, 对三块软骨分别进行区域的区域生长, 因此, 区域生长的结果中可能会存在软骨重叠等状况, 需要对区域生长结果进行后处理, 保证软骨分割的准确性。根据解剖学知识, 髌软骨右侧与股软骨左侧会出现重叠现象, 股软骨下侧与胫软骨上侧会出现重叠现象。对于重叠部分, 分别计算该区域内体素点与相应两条邻近边缘的欧氏距离, 将重叠像素点判定为距离小的一类软骨。然后, 对软骨结果采用形态学的闭操作、孔洞填充等操作, 得到最终的软骨分类结果。

3 实验结果与分析

3.1 评价参数

本文给出了金标准及自动分割结果的二维切片形式, 可以对它们之间的差异进行视觉评价, 此外, 除了视觉评价, 还使用灵敏度(Sensitivity)、特异度(Specificity)、Dice相似度系数(DSC)来定量地比较本文自动软骨分割结果与金标准的差异。使用 $\begin{matrix} TP 、 TN 、 FP 、 FN \end{matrix}$ 分别表示真阳性、真阴性、假阳性、假阴性, 则以上3个评价标准的表达式为:

$\begin{matrix} \begin{matrix} Sensitivity = \frac{TP}{TP + FN} (9) \\ Sepicficity = \frac{TN}{FP + TN} (10) \\ DSC = \frac{2 \times TP}{2 TP + FP + FN} (11) \end{matrix} \end{matrix}$

式中:DSC为评价的重要指标, 代表金标准与自动软骨分割结果的空间重叠指数, 范围在0~1之间, 1代表完全重叠, 0代表无重叠。

3.2 结果与分析

手工分割的金标准、文献[14]方法与本文方法的软骨分割结果对比如图5所示。从图中可以看出, 与文献[14]相比, 由于优化了初始骨-软骨边缘, 本文优化了股软骨分割结果, 去掉了股骨左上方多余的软骨体素, 以及胫骨左右两端的软骨体素。而第7层切片由于采用文献[14]的特征FV2未能识别出髌软骨, 因此造成了软骨结果的缺失。与文献[14]相比, 本文方法在第15层面优化了髌软骨, 平滑了胫软骨。

	Figure Option View Download New Window
	图5 分割结果比较Fig.5 Result comparison the with different methods

采用本文方法及文献[14]方法对5个人的128幅图像进行定量分析, 结果如表3所示。从表3中可以看出, 与文献[14]方法相比, 本文方法有明显提高。例如, 第7层中, 文献[14]未分割出的髌软骨被本文方法准确分割出来。通过本文方法, 股软骨、胫软骨、髌软骨的平均DSC系数均达到0.8以上。从所呈现的视觉和定量分析结果表明, 本文方法能够提供有效的软骨分割结果, 且对于软骨分割精度有显著提升。

表3 使用不同方法的平均DSC、平均特异度和平均敏感度比较 Table 3 Mean DSC, mean sensitivity, mean specificity with different methods

4 结束语

提出了一种基于多特征SVM边缘定位和弹性区域生长的膝关节软骨自动分割算法。首先, 通过改进的自适应Canny算子获取初始边缘。然后, 提取初始边缘的多种特征, 并采用一对一的多分类SVM算法进行分类, 实现骨-软骨边缘的定位。最后, 在定位的骨-软骨边缘附近自适应地选取种子点, 并采用弹性区域生长方法进行阈值分割, 从而完成对膝关节软骨的准确分割。本算法的视觉和定量分析均表明其分割结果较相关文献明显提高, 且与手工分割金标准相近, 显示其能够实现软骨的精确分割, 为后续三维建模和计算形态参数等提供了数据基础。

下一步将对提取的边缘特征进行筛选, 从而实现更精确的边缘定位及最终软骨分割结果。此外, 对于多层MRI序列, 相邻层之间AC所占面积相近、对比度相近、在图像中的位置相近, 可以考虑利用层间先验知识来为下一层分割提供先验区域, 将可能明显减少时间代价, 并排除先验区域之外的干扰, 从而提高分割精度和稳定性。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Hani A F M, Kumar D, Malik A S, et al. Multinuclear MR and multilevel data processing: an insight into morphologic assessment of in vivo nnee articular cartilage[J]. Academic Radiology, 2015, 22(1): 93-104. [本文引用:1]
[2]	Peterfy C G, Guermazi A, Zaim S, et al. Whole-organ magnetic resonance imaging score (WORMS) of the knee in osteoarthritis[J]. Osteoarthritis & Cartilage, 2004, 12(3): 177-190. [本文引用:1]
[3]	Hani A F M, Kumar D, Malik A S, et al. Non-invasive and in vivo assessment of osteoarthritic articular cartilage: a review on MRI investigations[J]. Rheumatology International, 2015, 35(1): 1-16. [本文引用:1]
[4]	Hani A F M. Fusion of multinuclear magnetic resonance images of knee for the assessment of articular cartilage[C]∥International Conference of the IEEE Engineering in Medicine and Biology Society, Osaka, 2013: 6466-6469. [本文引用:1]
[5]	Jaremko J L, Cheng R W, Lambert R G, et al. Reliability of an efficient MRI-based method for estimation of knee cartilage volume using surface registration[J]. Osteoarthritis & Cartilage, 2006, 14(9): 914-922. [本文引用:1]
[6]	König L, Groher M, Keil A, et al. Semi-automatic segmentation of the patellar cartilage in MRI[J]. Informatik Aktuell, 2007, 25(27): 404-408. [本文引用:1]
[7]	Swanson M S, Prescott J W, Best T M, et al. Semi-automated segmentation to assess the lateral meniscus in normal and osteoarthritic knees[J]. Osteoarthritis & Cartilage, 2010, 18(3): 344-353. [本文引用:1]
[8]	Shim H, Chang S, Tao C, et al. Knee cartilage: efficient and reproducible segmentation on high-spatial-resolution MR images with the semiautomated graph-cut algorithm method[J]. Radiology, 2009, 251(2): 548-556. [本文引用:1]
[9]	Lee J G, Gumus S, Moon C H, et al. Fully automated segmentation of cartilage from the MR images of knee using a multi-atlas and local structural analysis method[J]. Medical Physics, 2014, 41(9): 092303. [本文引用:1]
[10]	González G, Escalante-Ramírez B. Knee cartilage segmentation using active shape models and local binary patterns[J]. Optics Photonics & Digital Technologies for Multimedia Applications III, 2014, 9138: 91380K. [本文引用:1]
[11]	Parsoon A, Petersen K, Lgel C, et al. Deep feature learning for knee cartilage segmentation using a triplanar convolutional neural network[J]. Medical Image Computing and Computer Assisted Intervention (MICCAI), 2013, 8150: 246-253. [本文引用:1]
[12]	Parsoon A, Petersen K, Lgel C, et al. Femoral cartilage segmentation in knee MRI scans using two stage voxel classification[C]∥Conf Proc IEEE Eng Med Biol Soc, Osaka, 2013: 5469-5472. [本文引用:1]
[13]	庞剑飞, 邱明国, 陈伟, 等. 基于边缘检测与支持向量机的关节软骨自动分割算法研究[J]. 第三军医大学学报, 2013, 35(16): 1653-1657. Pang Jian-fei, Qiu Ming-guo, Chen Wei, et al. Segmenting articular cartilage automatically by edge detection and support vector machine[J]. J Third Mil Med Univ, 2013, 35(16): 1653-1657. [本文引用:1]
[14]	Pang J, Li P Y, Qiu M, et al. Automatic articular cartilage segmentation based on pattern recognition from knee MRI images[J]. Journal of Digital Imaging, 2015, 28(6): 695-703. [本文引用:7]
[15]	Liu H, Jezek K C. Automated extraction of coastline from satellite imagery by integrating Canny edge detection and locally adaptive thresholding methods[J]. International Journal of Remote Sensing, 2010, 25(5): 937-958. [本文引用:1]
[16]	Gonzalez R C, Woods R E. Digital Image Processing[M]. 3rd ed. US: Prentice Hall, 2009: 784-788. [本文引用:1]
[17]	Theodoridis S, Koutroumbas K. Pattern Recognition[M]. 4th ed. Greece: Elsevier, 2009: 134-137. [本文引用:1]
[18]	Hsu C W, Lin C J. A comparison of methods for multi-class support vector machines[J]. IEEE Transactions on Neural Networks, 2002, 13(2): 415-425. [本文引用:1]
[19]	Hojjatoleslami S A, Kittler J. Region growing: a new approach[J]. IEEE Transactions on Image Processing, 1998, 7(7): 1079-1084. [本文引用:1]

2015

0.0

... 研究表明^[1],OA的发作和退化等级与膝关节软骨(AC)密切相关 ...

2004

0.0

... 由于AC的损害程度反映疾病的严重程度,因此临床常用AC的损害程度作为诊断和评价OA退化程度的重要指标^[2] ...

2015

0.0

... 磁共振图像(MRI)具有无创性、无辐射性以及良好组织分辨率等性质,已经被证明是监测AC损害程度(即形态及生理变化)的最适宜的成像检测手段^[3,4] ...

2013

0.0

... 磁共振图像(MRI)具有无创性、无辐射性以及良好组织分辨率等性质,已经被证明是监测AC损害程度(即形态及生理变化)的最适宜的成像检测手段^[3,4] ...

2006

0.0

... 目前,临床上通常采用全手动或半自动的AC分割方法^[5,6,7,8],这些方法费时费力、效率较低,因此设计一种面向AC精准、快速的全自动分割方法显得尤为迫切 ...

2007

0.0

2010

0.0

2009

0.0

2014

0.0

... Lee等^[9]使用图谱和局部结构分析方法提取关节软骨 ...

2014

0.0

... González等^[10]提出了一种基于主动轮廓模型(ASM)与局部二进制模式(LBP)相结合的半自动分割方法 ...

2013

0.0

... Prasoon等^[11,12]提出了3个二维卷积神经网络(CNNs)组成的新型像素分类整合系统,这3个网络分别对应三维图像的 xy、yz和 zx平面 ...

2013

0.0

... Prasoon等^[11,12]提出了3个二维卷积神经网络(CNNs)组成的新型像素分类整合系统,这3个网络分别对应三维图像的 xy、yz和 zx平面 ...

2013

0.0

. 2013, 35(16):1653-1657

Segmenting articular cartilage automatically by edge detection and support vector machine

基于边缘检测与支持向量机的关节软骨自动分割算法研究

Pang Jian-fei , Qiu Ming-guo , Chen Wei

庞剑飞, 邱明国, 陈伟

目的设计并实现一种针对核磁共振成像(MRI)图像的关节软骨自动分割算法。方法利用像素的整体与局部特征分别构建二分类支持向量机(support vector machine,SVM)分类器对股软骨、胫软骨及髌软骨进行自动分割。首先提出一种基于边缘数目反馈的Canny检测器阈值迭代法并利用该方法提取图像的主要边缘,随后根据特征参数对提取的边缘进行识别并标记出不同的骨-软骨边缘,利用训练的SVM分类器对软骨进行初步分割并根据软骨的解剖位置缩小搜索空间,最后利用形态学操作对初步分割结果进行优化。结果自动分割结果中软骨的形态轮廓与原始图像吻合效果好,股软骨、胫软骨及髌软骨的Dice’s系数平均值分别为0.80、0.76、0.74,与手工分割结果具有较好的一致性。结论该算法能够准确、快速地分割出MRI图像中不同的软骨组织。

... 庞剑飞等^[13,14]使用边缘定位与支持向量机(SVM)相结合的方法自动分割软骨 ...

2015

0.0

... 庞剑飞等^[13,14]使用边缘定位与支持向量机(SVM)相结合的方法自动分割软骨 ...

... 2 结果与分析手工分割的金标准、文献[14]方法与本文方法的软骨分割结果对比如图5所示 ...

... 从图中可以看出,与文献[14]相比,由于优化了初始骨-软骨边缘,本文优化了股软骨分割结果,去掉了股骨左上方多余的软骨体素,以及胫骨左右两端的软骨体素 ...

... 与文献[14]相比,本文方法在第15层面优化了髌软骨,平滑了胫软骨 ...

... 采用本文方法及文献[14]方法对5个人的128幅图像进行定量分析,结果如表3所示 ...

... 从表3中可以看出,与文献[14]方法相比,本文方法有明显提高 ...

... 例如,第7层中,文献[14]未分割出的髌软骨被本文方法准确分割出来 ...

2010

0.0

... 但是,传统Canny算子的高、低阈值 Th和 Tl需手动设置,而膝关节MRI图像的切片层数在25层左右,手动设置 Th和 Tl不仅耗时,而且单纯地以人的感官选取会存在主观性强、阈值普适性差等问题,因此本研究采用自适应Canny边缘检测算法^[15],动态选取合适的 Th和 Tl以提高后续的边缘检测效率和普适性 ...

2009

0.0

... Step2 对 ImgradNMS采用迭代算法^[16]得到全局阈值 Tgrad,并计算 Tgrad左右两个区域的均值 μl和 μh,以及方差 σl和 σh,动态地设定高、低阈值为: ...

2009

0.0

... SVM^[17]是一种监督机器学习方法,将输入空间线性不可分的样本通过核函数映射到一个高维特征空间,转换为在特征空间中求解线性约束下的凸二次规划问题,目标是找到一个最优的分离超平面,将所有的训练向量正确分类 ...

2002

0.0

... 本研究采用一对一方法^[18]来解决该多分类问题 ...

1998

0.0

... 4 结合骨-软骨接触面的区域生长区域生长法^[19]的基本思想是:选取一组种子点,将与种子预先定义的性质相似的那些邻域像素添加到每个种子上来形成这些生长区域 ...