第 25 卷 第 6 期 2008 年 6 月 计 算 机 应 用 研 究 Application Research of Computers Vol. 25 No. 6 Jun. 2008 一 种 测 地 线 活 动 轮 廓 模 型 的 快 速 算 法 * 张 博, 苏永利, 张书玲 ( 西 北大 学 数学 系, 西 安 710069) 摘 要: 为 了完 成测 地线 活动 轮廓 模型 在图 像分割 中的 平滑 预处 理和 图像 梯度 的计 算, 给出 水平 集 方 法以 实 现 中符 号距 离函 数的构 造。 通过 对 Gaussian 函 数以 及差 分运 算 的讨 论 , 提 出 了 一种 基 于 Sobel 算 子图 像 预 处 理 方 法, 并利 用同 心圆 扩散方 法计 算符 号距 离。 得出 将图像 平滑 与梯 度计 算结 合为 Sobel 算子 的处 理, 一种 计算 符 号 距离 函数的 快速 计算 方法 。通 过实 验比较 , Sobel 算子既 可以 完成 平滑 处理 , 又 降 低 了差 分 计 算梯 度 的 时间 复 杂 度, 而同 心圆 扩散 方法 能够 提高 模型算 法的 执行 效率 。 关键 词: 水 平集 ; 测地 线活 动轮 廓; Gaussian 函 数; Sobel 算 子; 符 号距 离函 数 中图 分类 号: TP301 文 章编 号: 1001- 3695( 2008) 06- 1765- 03 文 献标 志码: A Fast algorithm of geodesic active contour ZHANG Bo, SU Yong-li, ZHANG Shu-ling ( Dept. of Mathematics, Northwest University, Xi’an 710069, China) Abstract: To smooth the image and calculate the gradient of image in the process of using GAC Model in image segmentation, and to construct the signed distance function in level set method. By discussing the calculation process of Gaussian function and difference approximation, a new pre-process method of GAC based on Sobel operator was presented. Then by expending of circle, the paper presented a method to calculate the signed distance. The operator proposed could smooth the image and cal- culate the gradient simultaneously and the method could calculate the distance quickly. The results of experimental comparison show that the algorithm proposed both enhanced smooth effect and reduced the time complexity in completing the GAC. Key words: level set; geodesic active contour( GAC) ; Guassian function; Sobel operator; signed distance function( SDF) 0 引言 G . Sapiro[ 1] 提出的测地线活动轮 廓模型 ( GAC) 是目前 基 于偏微分方程( PDE) 方 法图像 分割的 最基本 的活动 轮廓模 型 之一。该模型将图像分割问题 归结为 极小化 一个封 闭曲线 的 能量泛函, 并利用变分法将极小化能量泛函转换为关于封闭曲 线的梯度下降流; 然后利用 PDE 方法 完成曲线演 化, 并使 演化 过程在对象边缘处停止, 从而完成对图像的分割。为提高该模 型的分割效果以及算法的执行速度, 本文主要进行了如下两方 面的工作: a) 由于该方法在实现过 程中, 需要 计算图 像的梯 度模, 而 实际的图像往往存在噪声或伪边缘, 会使得计算出的梯度模比 较杂乱, 对后续的曲线演化结果影响较大。有必要在计算之前 对图像进行一定的平滑预处 理, 以消除 噪声或 伪边缘 的影响。 本文采用 Gaussian 函数进行图像的平滑处理, 然后再利用中心 到闭合曲线( 水平集曲线) 的符号距 离, 如果 直接计算, 计算 复 杂度为 O( MN) [ 1] 。其中: M 为 C 上 总点 数; N 为图 像的 总 像 素数, 可见计算量很大。Malladi 提出的快速步进 法( fastmarch- ing method) 生成的 SDF[ 1] , 计算量为 O( N ln M) 。本文首先 对 每个像素点邻近区域曲线形状进行分析, 快速确定曲线的内外 点; 再对图像像素点之间的 距离特 征进行 分析, 并构 造距离 表 ( 该距离表在以后 曲线 演化 的每 一步 不再改 变) ; 然后 对每 一 个像素点通过以该点为圆心的同心圆不断向外扩散, 找到与曲 线相切的同心圆, 该同心圆的半径就是圆心像素点到曲线的距 离; 最后查找距离表, 即得 圆心点 到曲线 的距离。该方 法计 算 像素点到曲线距离的计算复杂度为 O( N) , 而且 此算法的计 算 复杂度与闭合曲线所含点数无关。 1 测地线活动轮廓模型 G . Sapiro 在文献[ 1] 中提出了如下的能量泛函:  差分计算梯度, 并 将 两者 结合 为 Sobel 算 子, 利 用该 算子 可 同 LR = ∫L( C) 0 g( | I |[ C ( s ) ] ) ds ( 1 ) 时进行平滑预处理和梯度计算。 其中: s 为曲线的弧长参数; L( C) 表示曲线 C( s) 的总长; I 为 待 b) 在 GAC 模型的水平 集方法 实现中, 为保 持水平 集函 数 分割的图像; 函数 g 可选择 为任一 单调递 减函数, 于 是极小 化 ( 通常选择符号距 离函 数 SDF) 的特 性, 在每 次 迭代 后需 要 对 式( 1) 可理解为寻找 接近 于对 象边 缘并 以 g 为加 权函 数的 最 SDF 进行初始化。初始化方 法的实 现效率 将对整 个模 型的 实 短程闭合曲线 C。利用变分法极小化式( 1) 可得其对应的梯度 现过程产生重要的影响。SDF 的构 造需要 计算图 像各 像素 点 下降流为 收 稿日期 : 2007- 06- 15; 修 回日期 : 2007- 09- 05 作 者简介 : 张博( 1971- ) , 男, 陕 西西安 人, 博士 , 主 要研究 方向 为图像 处理 ( bobozhang000@ 163. com) ; 苏永利 ( 1974 - ) , 男 , 陕 西咸 阳人, 博士 , 基 金项 目: 陕 西省 教育厅 专项 资助项 目( JK05303) 主要研 究方向 为图 像处理 ; 张书玲 ( 1951 - ) , 男, 教授 , 博 导, 主要 研究方 向为 图像处 理. 

·6671· 计 算 机 应 用 研 究 第 25 卷 C / t = gkN - (  g ×N ) N ( 2) 但是, 速度场的自然延拓却导致函数 u 在经过演化后不再成为 其中: g = g( | I|) ; k 为曲线的曲率; N 为曲线 C 的内法向量。 距离函数。为了保持函数 u 作为距离函数的特性, 在每次迭代 式( 2) 的水平集方 法实现 中, 可引入 嵌入函 数 u = u( x, y, 后, 需要对距离函数重新初始化。 t) , u0 = u( x, y, 0) = C0, 通常将 u 选 做曲线 C 对应符 号距离 函 数, 即 SDF 的建立分为以下两部分: a) 首先建立符号图, sign P = ±1, P 为图像中任一点, 用于 区分 闭合曲 线的 内外 部, 并以 正 d [ ( x , y ) , C ] ( x, y) 在 C 内 号( 负号) 标记曲线 的内部 ( 外 部) ; b) 计算图 像中每 一像素 点  { u( x, y) = 0 从而式( 2) 对应的偏微分方程为 [ 3] - d [ , C ] y ) x , (   u/ t = | u |div g u/| u |   其中: div( ·) 为向量的散度。 2 预处理 ( x, y) 在 C 上 ( x, y) 在 C 外 P( x, y) 到曲线 C 的最短距离 d( P, C) , 则 u( x, y) = sign P ×d( P, C) ( 7 ) 3. 1 曲线内外点的判断 ( 3) 目前, 有关符号距离函 数构造 的文章, 对曲 线的内 外点 的 判断方法很少提到, 但内外点的判断也是影响曲线演化效率的 在利用式( 3) 进行曲线演化的 过程中, 需要计算 原始图 像 的梯度 I, 但是如果直接利 用原 始图像 进行 计算 所得 的模 值 图往往十分杂乱, 从而使 得曲线 在演化 过程中 不能准 确定位, 因此有必要在进行差分计 算梯度 前对图 像进行 平滑处 理。为 此, 笔者选择 Gaussian 函数对图像进行滤波: Iσ( x, y) = I( x, y) ×gσ( x, y) ( 4) 其中: gσ( x, y) = ( 1/σ 2πexp [ 槡 - ( x2 + y2) /2σ2] 。 为计算 方便, 通 常选 择半 窗宽 σ= 1 的 Gaussian 函 数, 从 而 Iσ( x, y) 图像在点( i, j) 处关于 变量 x 的一 阶差 分可 如下 计 算:     Iσ/ x ( i, j) = ( 1 /3 ) [ ( I/ x   I/ x ( i,   ( i, j - 1) + 2( I / x) j+1 ) ] ) ( i, j) + ( 5) 图像 I 的导数可采用中心差分计算, 即   I/ x ( i,   I/ x   I/ x ( i, j- 1) = [ I ( i + 1, j - 1 ) - I( i - 1 , j - 1) ] /2 j- 1) = [ ( ( i, j+1) = [ I I ( i + 1, j) - I( i - 1 , j ) ] /2 i + 1, j + 1 ) - I( i - 1 , j + 1) ] /2 一个关键因素。 任何一个闭合曲线不论其形状变化有多复杂, 在某一像素 点的附近其形状可分为四类( 图 3) 。 上下型 下上型 下下型 上上型 图3 像素点附近曲线形状 当人们从图像中某一 行的第 一个像 素点( 可 以认 为该 点 在曲线之外) 沿该行 从左 向右 进 行扫 描时, 如果 遇到 上 下型、 下上型两类边界时, 在该点 处发生 内外点 的转换; 如 果遇到 下 下型、上上型两类边界时, 则不发生转换。于是, 当进行逐行 扫 描时, 可通过判断遇到的边界的类型, 对该行的像素点作标记, 从而给出曲线内外点的判断, 具体实现步骤如下: a) 对图像作二值化处理, 使边界 点的像 素值为 0, 其余 点 像素值为 1。 b) 标记每一行的首元素为 out. ( 即 外点. 必要时 可将图 像 向左延拓一列, 再标记) 。 c) 逐行扫描。 代入式( 5) 并化简为   Iσ/ x ( i, [ I j) = ( 1 /6) I( i + 1, j + 1) ( i + 1, j - 1 ) + 2 ×I( i +1 , j) + - I( i - 1, j - 1) - ( a) 若 u( i, j) = 0( 边 界点 ) , 分 别判 断其 八邻 域点 中上 面 三个点或下面三个点中是否有边节点, 并用 up、down 记录。 若上面三个点中有边界点, 令 up = 1, 否则 up = 0; 2 ×I( i - 1, j) - I( i - 1, j + 1) ] ( 6) 若下面三个点中有边界点, 令 down = 1, 否则 down = 0; 式( 6) 的矩阵 形式 如 图 1 所 示, 即 权 系数 为 1 /6 的 Sobel 算子, 用此算子作用 于图 像 I, 即 得 Iσ 关于 变 量 x 的 差 分图。 同理, Iσ 关于变量 y 的差分图也可表示为图 2。 ( b) 若 u( i, j) = 1( 非 边界点) , 判断 该点 的前 一个 像素 点 ( i, j - 1) 是否为非边界点: ①若非边界点, 则该点 与其前 一像素 点做相 同标记, 进 行 -1/6 -1/3 -1/6 0 0 0 1/6 1/3 1/6 图 员 关于 x 的Sobel 算子 -1/6 -1/3 -1/6 0 0 1/6 1/6 0 1/3 图2 关于 y 的 Sobel 算子 3 符号距离函数的建立 GAC 模型的水 平集 方 法 实现 中, 式 ( 2) 是 针 对曲 线 的 演 化, 故其演化速度 V = gk -  g ×N 只 在曲 线 C 上 有意 义。而 当由式( 2) 得到水 平 集函 数的 演 化式 ( 3) 时, 演化 速度 Vext = div( g u| u|) 却需要在整个图像上有定义, 这就需要对 速度  下一像素点的判断; ②若是边界点, 判断 up ×down 的值: 若 up ×down = 1, 则 遇到图 3 中 前两 种 边界, 此 时该 点 标 记与其所在行中左侧距离该点最近的非边界点标记相同; 若 up ×down = 0, 则 遇到图 3 中 后两 种 边界, 此 时该 点 标 记与其所在行中左侧距离该点最近的非边界点标记相异。 ( c) 令 up = 0, down = 0, 进行下一像素点的判断。 图 4 是利用以上方法生成的符号图, 整个算法只需对图 像 作一次扫描即可完成。 3. 2 距离函数的建立 函数重新定义, 通常使用的方法是速度场的自然延拓, 即 图像上任意像素点 P0( x0, y0) 到曲线 C 的距离为 Vext ( x, y) = V( x, y) ( x, y) ∈C d( P0 , C) = min { d( P0 , P) |P∈C } 

第 6 期 张 博 , 等 : 一种 测地线 活动 轮廓 模型 的快 速算 法 ·7671· 3. 2. 1 距离函数表的建立 计算精度高, 且计算时间随着曲线点数的增加会缩短。 对于任一网格 点, 它 到其 他网 格 点的 距离 为 一系 列 离 散 值。如图 5 所示, 以网格点 P0 为 圆心的 一系列 同心圆 的半 径 就是点 P0 到圆上各网格点的距离; 如图 6 的矩阵所示, 矩阵中 只给出了 1/8 邻域中网 格点的 距离值, 其余 网格点 到 P0 的 距 离可利用对称关系得到。 (a)初始曲线 (b)符号图 图 4 初始曲线及其符号图 02+12 员2+12 0 2 5 2 1 02+22 员2+22 圆2+22 02+32 噎 员2+32 噎 圆2+32 噎 猿2+32 噎 左 噎 图 远 网格点间距离表渊dist) 图5 同心圆扩散 3. 2. 2 同心圆扩散 对网格点 P0, 每次判断以 P0 为圆心的每 个同心圆上 的网 格点中是否有曲线上的点; 若没 有则考 察下一 个更大 的圆, 直 至圆上的网格点有曲线上的点为止。此时, 该圆的半径就是点 P0 到曲线的距离, 直接在距离表中查出该距离值即可。 具体编程时, 只 需要 分析 点 P0 的 1 /8 邻 域, 邻 域中 其 余 的点与 1 /8 该 邻 域 中 的 点 分 别 关 于 x 轴、y 轴、坐 标 原 点、 y = ±x对称, 即若以 P0 为圆心的同心圆上有点( i + s, j + k) , 则 一定有点( i ±s, j + k) 和( i ±k, j ±s) ; 为 避免编 程时, 邻域点 的 下标范围超出图像, 可将 原始图 像作一 定延拓, 而不 影响运 算 量, 具体算法如下: a) 将图像 u( 大小为 M ×N) 二值化, 边界点的像素值 为 0, 其余为 1。 b) 将图像 对称 延拓为 U( 大 小为 3M ×3N) , 使 U( M + 1: 2M, N + 1: 2N) = u; c) 遍历所 有 网 格 点, 计 算 每 一 个 网 格 点到 曲 线 的 距 离。 对网格点 P0( i, j) 作同心圆扩散: s = 1: M; k = 0: s。 网格大小 表员 128伊128 256伊256 缘员圆伊缘员圆 128伊128 256伊256 缘员圆伊缘员圆 128伊128 256伊256 缘员圆伊缘员圆 杂阅云 构造方法比较 282/839 505/1890 1390/2668 282/839 505/1890 1390/2668 282/839 员愿怨园 曲线点数/两种 计算用时/s 7.45/28.32 60.0/177.8 722/1186 源.27/15.11 36.5/90.2 397/605 圆.81/10.53 20.3/49.1 203/357 员猿怨园辕圆 6远愿 计算方法 直接法 快速步进法 同心圆扩散法 4 实验比较 图 8 和 9 是对同一幅图像分别没有经过 Sobel 算子预处理 和使用 Sobel 算子预 处理 后 的分 割结 果。从 分 割结 果可 以 看 出, 进行预处理后, 对对象轮廓的定位比较准确。 (a)0 (b)500 (c)2000 (d)4000 图 8 曲线演化(没有使用 Sobel 算子) (a)0 (b)500 (c)2000 (d)4000 图 9 曲线演化 (使用 Sobel 算子) 5 结束语 本文给出了一种 利 用 Sobel 算 子对 测地 线 活动 轮廓 模 型 进行预处理的方法。该方法很 好地消 除了噪 声或伪 边缘对 曲 线演化结果的影响。由于将预 处理和 计算梯 度两个 工作结 合 为一个算子, 使得处理过程容易实现, 提高了运算速度。同时, 借助于同心圆扩散方法构造符号距离函数, 也使得曲线的演化 速度大大提高, 从而 为基 于 PDE 的 图像 分割 方法 提供 了有 力 的支持。 若 U( i ±s, j ±k) 、U( i ±k, j ±s) 其中之一为 0, 记为点 P, 则 参考文献: d( P0 , C) = d( P0P) , 该距离对应距离表中的 dist( s, k) 。 d) 建立 SDF: d = signP0 ×d( P0, C) 。 图 7 是利用以上算法生成的符号距离函数示意图。 (b) (a) 图 7 原始图像与符号距离函数 [ 1] APIRO G. Geometric partial differential equations and image analysis [ M] . New York: Cambridge University Press, 2001. [ 2] 陆金甫 , 关治 . 偏 微 分 方 程 数 值 解 法 [ M] . 北 京: 清 华 大 学 出 版 社, 2004 . [ 3] SETHIAN J A. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vi- sion and materials science [ M] . Cambridge, United Kingdom: Cam- bridge University Press, 1999. [ 4] SETHIAN J A . Curvature and the evolution of fronts [ J] . Communi- cation of Mathematical Physics, 1985, 101( 4) : 185- 197. 表 1 对常用的三种 SDF 方法( 直接法、快速步进法、同心圆扩 [ 5] COHEN L, KIMMEL R. Globel minimum for active contour models: 散法) 进行了试验对比。从实验结果可以看出, 在不同网格大小 a minmal path approach[ J] . Internatio nal Jo urna l of Co mputer Vi- ( 128 ×128, 256 ×256, 512 ×512) 和不同的曲线长度下, 本文算法的 sio n, 1997, 24 ( 17 ) : 57- 78.