随着科技的进步, 医学影像技术成为当今诊断疑难杂症的重要手段, 其中, CT成像技术发展较早, 技术相对成熟, 在临床的各个领域获得了较为广泛的应用.其原理是利用人体不同组织对X光吸收率不同, 通过接收转换衰减信号获得人体各部分的医学图像^[1].但是人体腹腔环境复杂, 器官与主动脉密度相近, 目前没有一个理想的腹腔CT主动脉分割算法, 很多时候还得依靠人工分割, 造成工作强度大, 分割效果差, 无法完整分割主动脉结构, 错分、伪影、断层等现象时有发生.对后期主动脉三维重建以及主动脉疾病的诊断和治疗埋下了隐患.

1988年, Osher等人^[1–3]提出了描述曲线演化过程的水平集方法(LeveI Set Methods).其基本思想是将当前低维度的曲线问题转化为更高维空间中的水平集函数曲线的隐式解.为了获得准确的边缘检测结果, 许多人希望在水平集演化过程中加入附加信息, 例如区域信息、形状信息和相位信息等. Belaid等人^[4,5]提出了一种基于相位的水平集分割算法: 利用单演信号所产生的局部相位和测量局部图像和轮廓的正常运动方向之间的局部相对方向对图像进行分割, 但是需要时刻调整相关参数. Estellers等人^[6,7]提出一种利用图像几何形状的水平集分割方法, 其利用相邻水平集梯度的相关性调整目标轮廓, 但对于灰度强度不均匀的图像, 其效果不佳. 为了解决水平集方法计算效率较低的问题, 1993年Chopp^[8]提出了窄带法以提高求解效率; 1995年Adalsteinsson和Sethian^[9]把演化计算限制在了一个零水平集函数附近的窄带范围内, 从而大大缩短了曲线演化速度, 并提高了其计算精度; 2001年Chan和Vese提出了一种性能较好的全局优化分割模型, 但只能应用于所需分隔物体与背景差别较大图像, 应用范围受大很大限制; 2005年Li^[10]等人进一步改进其算法, 对浓度不均的图像, 构造区域可求拟合能量模型, 得到了较好的分割结果; 为了进一步提高处理速度, Bresson等人提出基于全局变分思想的对偶公式迭代算法, 避开求解偏微分方程, 明显提高了演化速度. 在前人工作基础上, 考虑局部特征的强度与方向对水平集演化加权, 提出一种基于局部边缘特征加权水平集演化的主动脉CT图像分割方法. 本方法将演化方向相邻区域和演化轮廓内外的梯度方向同时进行了测量, 并考虑演化轮廓相邻区域的平均边缘强度, 对曲线面积和曲线长度赋予不同的权重, 使得弱边缘对分割精度影响最小, 实验结果表明, 此方法能够得到精确的边缘检测结果.

1 局部加权水平集演化算法

基于活动轮廓的变分水平集算法一般需要加入图像信息来定义目标函数^[11], 本算法使用边缘信息作为图像特征, 用不断变化的主动脉轮廓驱动算法寻找边缘, 流程图如图1所示.

边缘强度检测公式如下:

$g \triangleq \frac{1}{{1 + {{\left| {\nabla {G_\sigma }*I} \right|}^2}}}$

(1)

其中, $g \in \left[ {0, 1} \right]$ , I表示域Ω上的图像, ${G_\sigma }$ 是具有标准偏差的高斯核函数, 对二者进行卷积运算.基于函数g, 定义水平集(LSF)的能量函数 $\varphi $ :

$E\left( \varphi \right) = R(\varphi ) + Length(\varphi ) + Area(\varphi )$

(2)

水平集函数通常在零水平集附近过于陡峭或平滑会导致误差, 因此引入 $R(\varphi )$ 距离正则化项使其保持形状; $Length(\varphi )$ 表示要检测的曲线长度, $Area(\varphi )$ 表示要检测的曲线面积, 为了使边缘能量最小, 根据公式(1)分别定义 $Length(\varphi )$ 和 $Area(\varphi )$ 项如下^[12]:

$Length\left\{ {\varphi = 0} \right\} = \mathop \smallint \limits_\Omega g\delta \left( \varphi \right)\left| {\nabla \varphi } \right|dx$

(3)

$Area\{ \varphi \geqslant 0\} = \mathop \smallint \limits_\Omega gH\left( \varphi \right)dx$

(4)

其中, $\delta \left( \varphi \right)$ 表示狄拉克δ函数, 用来计算g函数零水平集的线积分; $H\left( \varphi \right)$ 表示阶跃函数, 用来计算演化轮廓内的能量. $Length(\varphi )$ 在物体边界处取最小值, $Area(\varphi )$ 通过在平滑处加快演化速度, 在边界复杂处减慢演化速度帮助寻找期望的边界. 此算法在强边缘情况下效果较好, 但是在灰度不均匀等弱边缘情况下可能出现漏检或错检, 尤其是主动脉医学图像中, 原始图像显示效果不佳, 此算法的缺点被进一步放大^[13].

根据主动脉图像的特点, 对于弱边缘需要提高图像的边缘检测准确性, 提出一种加权函数, 根据边缘内部和外部相邻区域的图像特征, 为 $Length(\varphi )$ 和 $Area(\varphi )$ 项分配不同的优先级.设加权因子为ω, 根据局部边缘特征确定其重要程度, 带加权因子ω的公式如下:

$Length'\left\{ {\varphi = 0} \right\} = \mathop \smallint \limits_\Omega g\left( {1 - \omega \left( {\varphi , k} \right)} \right)\delta \left( \varphi \right)\left| {\nabla \varphi } \right|dx$

(5)

$Area'\{ \varphi \geqslant 0\} = \mathop \smallint \limits_\Omega g\omega \left( {\varphi , k} \right)H\left( \varphi \right)dx$

(6)

k是一个常数, 决定获取局部特征的区域大小. 公式 $\omega \left( {\varphi , k} \right)$ 由下式给出:

$\omega \left( {\varphi , k} \right) = I{\left( {\varphi , k} \right)^{\left( {1 - \gamma \left( {\varphi , k} \right)} \right)}}$

(7)

其中, $I \in \left[ {0, 1} \right]$ , 为轮廓曲线C上由k常数确定的平均边缘强度, γ∈[−1, 1], 为C上某点法线方向与其GVF场中k相邻域的内积. 轮廓曲线C的相邻域定义如下:

$\vartheta \left( {{\rm{\varphi }}, m} \right) = {\rm{\delta }}\left( {\rm{\varphi }} \right)\left| {\nabla \varphi } \right| + m\vec N$

(8)

其中, m的大小、正负与其和零水平集相对位置有关, 具体含义如图2所示, 其中虚线圆表示零水平集, 外圈等高圆m从整数1开始, 内圈等高圆m从整数–1开始.

由以上基础公式可以得出边缘平均强度公式如下:

$I\left( {\varphi , k} \right) = \frac{1}{{2k}}\mathop \sum \limits_{m = 1}^k \left[ {\mathop \smallint \limits_{\rm{\Omega }}^{} g\vartheta \left( {{\rm{\varphi }}, m} \right)dx + \mathop \smallint \limits_{\rm{\Omega }}^{} g\vartheta \left( {{\rm{\varphi }}, - m} \right)dx} \right]$

(9)

同样地, 计算轮廓曲线C的法向量与沿着k相邻轮廓的GVF场的平均内积γ, 如下:

$\gamma \left( {\varphi , k} \right) = \frac{1}{{2k}}\mathop \sum \limits_{m = 1}^k \left[ {\mathop \smallint \limits_{\rm{\Omega }}^{} \overrightarrow {N, } \vec V\vartheta \left( {{\rm{\varphi }}, m} \right)dx + \mathop \smallint \limits_{\rm{\Omega }}^{} \overrightarrow {N, } \vec V\vartheta \left( {{\rm{\varphi }}, - m} \right)dx} \right]$

(10)

其中, $\vec V$ 表示GVF场, 当 $\overrightarrow N $ , $\vec V$ 两者方向相近时, 最终结果γ的值趋近于1.

用带有加权因子ω的公式(5)、(6)分别替换能量函数 $E\left( \varphi \right)$ 中的(3)、(4), 即可用ω根据局部特征调节 $Length(\varphi )$ 和 $Area(\varphi )$ 项的优先级, 当零水平集位于非光滑边缘时, $Length(\varphi )$ 项成为主要的驱动项, 当零水平集位于光滑边缘时, $Area(\varphi )$ 项成为主要驱动项. 最后, 对能量函数 $E\left( \varphi \right)$ 求极小值可以转化成求Euler-Lagrange方程的解.

2 实验结果及分析

为了验证本文算法的正确性和有效性, 利用本文提出的算法对主动脉CT图像进行分割, 并与传统水平集演化算法^[14]作对比, 实验数据均来自医院实际拍摄的患者腹腔主动脉CT图像. 本文是在 Intel i5(R)CPU 2. 60 GHz, 4 GB 内存, Windows 7 的环境下, 应用 MATLAB所做的实验.

如图3所示, 为采集到的原始图像, 图4为采用传统水平集演化算法分割的主动脉轮廓, 可以看出, 其误将腹腔其他组织尤其是肾脏边缘误分割为主动脉, 分割效果较差; 图5采用基于局部边缘特征加权的水平集演化算法, 取得了较为理想的效果.为了具体衡量算法的优越性, 选用DSC(Dice Similarity Coefficient)方法计算各算法的分割精度^[15], 计算方法如下:

${{DSC}}\left( {{{S}}, {{G}}} \right) = \frac{{2\left| {S \cap G} \right|}}{{\left| S \right| + \left| G \right|}}$

(11)

其取值范围为[0, 1], 表示算法S与标准G分割图像的相似性, 一般说来, G为人工标记的主动脉边缘, 当DSC取1时, 表示算法与标准分割完全一致, 效果最好. 具体分割精度如表1所示本文算法在精度与稳定性两个方面均优于传统水平集演化算法^[16].

3 结论与展望

本文采用基于局部特征加权的水平集分割算法, 并将其应用于较难分割的腹腔主动脉问题中. 针对传统水平集演化算法在弱边缘处效果较差的劣势, 分别考虑边缘长度和面积两项, 通过加权因子w对不同边缘处采取不同的策略, 实现能够较好分割腹腔主动脉的目的. 实验结果表明, 本文算法在分割精度和稳定性方面比传统水平集演化算法更胜一筹, 为后续主动脉3D模型重构奠定了基础, 也为腹腔主动脉分割算法提供了新思路.