基于三维CT图像的骨表面检测算法

姚明　程远志

摘 要：为了更好地获取三维CT图像中的骨表面，本文利用轮廓跟踪的思想提出了一种三维骨表面检测算法。该算法首先通过三维轮廓跟踪方法找出初始的三维表面轮廓，然后利用一阶导数的方法对滤波去噪后的图像中的初始表面体素计算法线方向；为了保证该法线方向的足够精确，本文提出了一种预估-校正的方案，来对之前计算出的法线方向进行修正，从而得出一个更精确的法线方向；最后，利用此法线方向建立一维信号、进行边缘检测。实验结果表明本文算法有效的增强了轮廓跟踪算法的抗噪性，也较好的提取出来完整的骨表面。

关键词：医学图像处理；轮廓跟踪；法线方向校正；边缘检测

中图分类号：TP391.41 文献标识号：A 文章编码：2095-2163（2015）03-

An Algorithm for Bone Surface Detection in 3D CT Images

YAO Ming， CHENG Yuanzhi

（School of Computer Science and Technology， Harbin Institute of Technology， Harbin 150001， China）

Abstract： In order to get a better bone surface of 3D CT images， this paper proposes a 3D bone surface detection algorithm with the idea of contour tracing. First， this algorithm get an initial surface with 3D contour-tracing algorithm. Then， through the first derivative of filtered images， normal directions of initial surface voxels can be estimated. To insure the accuracy of the estimation， this paper proposes an estimation-correction scheme to correct the normal direction and get the corresponding better values. Finally， with the normal direction， a 1D signal can be constructed for edge detection. The experiment result shows that this algorithm enhances the anti-noisy ability of contour tracing algorithm and gets a better bone surface.

Keywords： Medical Image Processing； Contour Trace； Normal Direction Correction； Edge Detection

0 引 言

对于当代的医学来说，用计算机来处理各种医学图像是一项重要的自动化医疗技术。而在医学领域，分割CT图像则是引导手术中不可或缺的组成部分。相比传统的二维CT图像，三维CT图像可以通过旋转更加立体、直观地显示体骨的各部分特征，对临床诊断和治疗成效具有重大的裨益作用。由上可见，三维CT图像中部位的自动分割即已成为近几年计算机处理医学图像中的研究攻关重点，而由于骨CT图像还掺杂着其他组织，所以CT图像中骨表面的检测也是其全程实施和付诸实现过程中的一个难点。

对于二维图像，当今国内外已有许多较为成熟的二维图像的分割算法[1]。这些算法大部分都是基于分割物体为同质的研究假设，而后利用全局或者局部的阈值[2-4]的方法来进行分割。另外一种当下更为流行的方法就是动态轮廓模型的方法[5]。这种方法在医学图像处理中的实用性已然获得了确切证明[6-7]。此外，一般的边缘检测或者说是边缘过滤也是一套常见的CT图像骨分割方法。此种方法就是立基于在图像中边缘的像素强度变化要快于非边缘区域这一研究结论之上的。

作为二维图像的延伸，在三维图像中，上述二维图像中的方法很多都可以扩展到三维图像中。一般情况下，三维图像中的边缘检测研究都是基于每一个二维切片的表面跟踪[8-10]。研究中，在进行跟踪时，就需要给出表面体素所符合条件的定义，常见方法则是基于梯度的三维图像的边缘检测算法。这些算法都可以看作是二维算法的延伸。只是这些算法在一定程度上却都忽略了三维图像的一些信息，而且复杂度较高。

基于此，本文提出了一种基于法线方向的边缘检测方法。首先，利用传统的边缘检测方法，如轮廓跟踪等方法获取初始边缘体素，然后计算初始边缘点的法线方向，并且为了获取更精准的法线方向，本文又提出了预估-校正的方案；当得到法线方向后在法线方向上建立一维信号，实现了边缘检测。

1 三维表面跟踪算法

首先，研究要利用轮廓跟踪获取一个初始轮廓。在三维图像的跟踪算法中，最重要的就是引入表面体素的判定条件，使得在跟踪图像时可以直接判断出单层的表面体素。在这里，本文引入了过零点的判断规则，并凭此定义了表面体素的判断条件。这个方法类似于LoG（高斯-拉普拉斯算子），但是LoG算子容易在边缘检测中产生多层结果。为了保证最终结果只有一层，研究提出将LoG算子和SDGD（高斯梯度方向二阶导滤波）融合在一起，定义了离散空间中的表面判定条件。首先要进行如式（1）所示的滤波。

（1）

其中， 是高斯函数在梯度方向上的二阶导数，I即为原图像。

在得到了这个滤波结果后，就可以用式（2）中的过零点条件进行判定。符合此式中条件的体素，就属于初始表面体素。

（2）

2 预估-校正方案计算法线方向

2.1 法线方向预估

为了预估法线方向，本文首先对图像进行了高斯滤波，然后对滤波后的图像计算一阶倒数即可得到法线方向，预估方法较为简单，此处不再赘述。研究关键则在于如何校正法线方向，下面即从二维图像入手展开对本文方法的论述。

2.2 二维法线方向校正方案

假设 为初始边缘点i的法线方向，而 就表示两个相邻的初始边缘点i-1和i之间的夹角。这样，由于髋关节的骨表面是平滑并且连续的，所以 和 之间的差就应当非常小。在理想情况下，如果轮廓是一个圆，那么这个差则如式（3）所示：

（3）

但是，在实际的骨表面（本文实验中使用髋关节图像）中并不会出现真正的圆。为此可以通过引入一个随机分布在0附近的变量 来修改式（3），得到式（4）：

（4）

正如上述研究中提到的，正常情况下 应该是很接近0的。但是已有研究估算的法线方向一定有误差，而且，法线方向的误差越大， 的值就越大。此时，如果指定一个变量 ，那么就可在精确的和估算的法线方向值之间建立起一种关系，具体如式（5）所示：

（5）

其中， 表示被校正过的法线方向，为了校正结果，可考虑利用已经校正过的点i-1来估算当前点i的值，也就是如式（6）所示：

（6）

这样，就可以利用式（3）和（6）式来重写式（5），具体可以得到：

（7）

之后，则可将如下两个基础公式 和 代入式（7）中，进一步地可得到：

（8）

研究至此，在单片二维图像中的法线校正方案就已做到了完整呈现。综上可知，对于点i，就可以利用已经校正过的邻居点i-1来完成对其的校正。而且，式（8）中的 就相当于点i的另一种法线方向的预估值，这个值就是利用了其相邻点的法线方向计算得到的。归纳可得，式（8）其实就是利用了一个因子 而将两种法线估算方案实现了结合。而关于式（8）中因子 应该如何选取，将会在3.1中进行详述。

2.3 三维法线方向校正方案

该方案若要应用到三维图像中，关于各个体素之间法线方向关系的表示就是一个基础研究问题。在此，研究中选择使用了球坐标系来衡量整个三维空间。这样，对于三维空间中的法线方向，在球坐标系中就可以用天顶角和方位角两个角度来表示，具体如图1 （a）所示。

（a）两个法线方向 （b） 平面z=0中方位角的校正 （c） 平面x=0中天顶角的校正

（a）2 normal directions （b） Projection of azimuth angle （c） Projection of zenith angle

图1 球坐标系中的法线方向与投影

Fig.1 Normal direction and projection in spherical coordinate

图1 （a）中， 和 是两个相邻体素的法线方向，其中 是方位角， 是天顶角。通过这种表示方法，即可将三维角度转换到两个不同的二维空间中，而后分别利用上文中所得的方法对其进行校正，如图1（b）和图1（c）所示。

如图1（b）和图1（c）所示，在这两个平面中，将分别使用式（8）来对方位角和天顶角进行校正，由此就得到了完整的三维图像中法线方向的预估-校正方案。

3 参数估计与边缘检测

3.1 三维法线方向校正方案

从上文可以看出，目前需要制定取值方案的参数有两个，一个是高斯一阶导中高斯滤波的噪声标准差 ，另一个是高斯一阶导所得计算结果的预估可信度 。

显然，高斯一阶导计算的法线方向 的取值也是随着 的改变而变化的，确定了一个最佳 值，也就同时确定了最佳的 。所以，研究尝试对体素 找到 对 变化的稳定区间。那么如果使用一个选定的区间 中的不同的 来计算其法线方向，计算后的结果变化又小于一个选定的值 时，就可判定为估算结果是稳定的，并且设定 ，即通常情况下， 可以大致定义为：

（9）

式（9）中， ， ， 表示该区间是满足条件要求的 的取值区间。上述条件是：该区间中所有的 取值估算得到的点 的法线方向之间的差值都小于或者等于选定的 ，并且 也是满足这一条件的最大区间。也就是说，在 区间外的 值所求的法线方向的差异都大于 。在此种方法下，就可以找到生存周期最长的 值，并将作为 的最终取值。

此外，研究还需要确定每个点最佳的 取值以便之后的边缘检测算子使用。毫无疑问，最佳的 就在 这个区间上，考虑到研究使用的是最大梯度的边缘检测算子，所以应该选择会使得 能够取得最大值的 作为点 的最佳 取值，具体如式（10）所示。

（10）

3.2 构建一维信号进行边缘检测

依照上述过程得到了法线方向，同时获得了各个体素的最佳 取值后，即以体素v为中心沿其法线方向上的体素的像素值构建成一个一维信号 。基于此，就只要选用一种边缘检测算子来进行一维滤波即可。在本文的实验中，为了方便计算，选用的是高斯一阶导数滤波，其优越性也已经获得了证明[11]。

这样，就完成了算法的整体流程，通过最终的边缘检测算子得到了寻求的表面体素，由此而实现了表面检测。

4 实验结果与结论分析

理论上，三维图像其实就是一组二维图像排列而成，因此大部分的三维分割算法都是由二维图像变换改进而来。相应地，本文的实验结果也即分为单片二维图像的结果和完整三维图像结果两部分来进行探讨与分析。需要说明的是，本文的实验都是基于髋关节的CT图像的。

4.1 二维图像

图2的左右对比展示了算法在二维切片上的效果。显然，算法在二维切片上分割出完整闭合的轮廓。

（a）原图 （b） 分割图

（a） Oringinal Image （b） Segementation image

图2 实验结果对比

Fig.2 Comparison of experimental results

图3（a） 和图3（b） 是本算法结果和单纯轮廓跟踪算法结果的局部对比。可以看到图3（b）中本文通过后续的法线校正和再次进行检测之后的结果对这种情况有了明显的改善，这就说明本文研究提高了轮廓跟踪算法的抗噪能力。而图3（c）和图3（d）则是在本算法与Canny边缘检测算子的局部结果的对比呈现。很明显可以看出，Canny算子的检测结果包含了除骨骼以外的其他组织的边缘，并且骨骼边缘还不一定是完整且闭合的，而本文的算法只找到了骨的边缘，滤除了其他组织，形成了完整闭合轮廓。这也就从根本体现了本文算法的优势。

（a） 轮廓跟踪初始结果 （b） 完整算法结果 （c） Canny算子结果 （d） 本算法结果

（a） Result of contour tracing （b） Result of the proposed algorithm （c） Result of Canny （d） Result of the proposed algorithm

图3 局部结果对比

Fig.3 Local results comparison

4.2 三维图像

图4（a）是使用Volview中包含的常见基本表面检测算法对髋关节CT图像处理后的结果，图4（b）是本研究算法处理的初步结果。

（a） Volview所含基本骨表面检测算法 （b） 本文处理方法

（a） Bone detection result of Volview （b） Result of the proposed algorithm

图4 三维髋关节骨表面检测结果对比

Fig.4 Result comparison of 3D hip joint surface detection

对比以上结果可以看出，由于本文的算法是基于跟踪的，所以只检测到了研究需要的髋关节，而不会和传统的算法一样把脊椎也显示包含于结果之中。此外，本文算法又将两侧的其他组织边缘也实现了过滤。

除了从图像上的直观比较外，还需要将结果与手动分割得到的黄金标准做一比较。研究分别测量了四组三维图像分割结果与黄金标准的Dice系数和平均曲面距离（Average Symmetric Surface Distance），并且和Canny算子与轮廓跟踪算法进行了结果比较，具体如表1所示。

表1 三种方法在三维髋关节CT图像上的精度评价

Tab.1 Accuracy evaluation of the three segmentation methods using 3D hip joint CT image

Data Set 轮廓跟踪 Canny算子 本文算法

Dice（%） ASD（mm） Dice（%） ASD（mm） Dice（%） ASD（mm）

Set 1 83.79 2.15 89.84 1.73 91.07 1.22

Set 2 85.46 1.93 90.03 1.49 90.55 1.47

Set 3 79.67 2.53 85.27 1.95 87.06 1.80

Set 4 81.24 2.37 86.30 1.84 86.77 1.82

从表1的结果可以看出，本文算法的结果跟黄金标准的Dice系数要高于其他两种算法，而ASD距离基本上则要小于其他两种方法。这就说明本文的算法对其余两种算法的提升在效果上更为明显，并且又剔除了其他组织，因而具有更为可观的综合算法效能。

5 结束语

本文提出了一种基于法线方向的边缘检测方法。先利用传统的轮廓跟踪方法，获取初始边缘点，然后计算初始边缘点的法线方向，并且利用预估-校正的方案校正法线方向；得到法线方向后在法线方向上建立一维信号，进行边缘检测。实验结果表明本文的算法可以获取完整闭合的三维骨表面，有效地提高了轮廓跟踪的抗噪性，并且对比Canny算子来说可以单独提取出骨的轮廓而滤去了其他组织。除此之外，对于轮廓跟踪和Canny算子来说，本算法也获取了较高的Dice系数和更小的ASD。

参考文献：

[1] PAL N R， PAL S K. A review on image segmentation techniques[J]. Pattern recognition， 1993， 26（9）： 1277-1294.

[2] KOHLER R. A segmentation system based on thresholding[J]. Computer Graphics and Image Processing， 1981， 15（4）： 319-338.

[3] REDDI S S， RUDIN S F， KESHAVAN H R. An optimal multiple threshold scheme for image segmentation[J]. Systems， Man and Cybernetics， IEEE Transactions on， 1984 （4）： 661-665.

[4] OTSU N. A threshold selection method from gray-level histograms[J]. Automatica， 1975， 11（285-296）： 23-27.

[5] LOBREGT S， VIERGEVER M A. A discrete dynamic contour model[J]. Medical Imaging， IEEE Transactions on， 1995， 14（1）： 12-24.

[6] MIGNOTTE M， MEUNIER J. A multiscale optimization approach for the dynamic contour-based boundary detection issue[J]. Computerized Medical Imaging and Graphics， 2001， 25（3）： 265-275.

[7] 吉宏伟， 何江萍， 杨新. 基于层次上下文活动轮廓的三维 CT 肝脏图像分割[J]. 生物医学工程学杂志， 2014（2）： 032.

[8] BHATTACHARYA P， WILD D. A new edge detector for gray volumetric data[J]. Computers in biology and medicine， 1996， 26（4）： 315-328.

[9] WANG P， CHENG L， LIN B， et al. Boundary surface detection in 3D biomedical images based on local gradient thresholds[C]//Biomedical Engineering and Informatics （BMEI）， 2011 4th International Conference on.， Shanghai： IEEE， 2011， 1： 175-179.

[10] PUDNEY C， ROBINS M， ROBBINS B， et al. Surface detection in 3 D confocal microscope images via local energy and ridge tracing[J]. Journal of Computer-Assisted Microscopy， 1996， 8（1）： 5-20.

[11] PAPLINSKI A P. Directional filtering in edge detection[J]. Image Processing， IEEE Transactions on， 1998， 7（4）： 611-615.