CN103065333A - 一种基于四元数的彩色snake图像分割方法 - Google Patents

Abstract

本发明属于计算机图像处理领域，为一种基于四元数的彩色snake图像分割方法。①输入一幅待分割的RGB空间彩色图像，并将该图像表示到四元数空间；②计算①中得到的四元数空间下表示的彩色图像的梯度幅值；③在待分割的RGB空间彩色图像上勾勒出感兴趣物体的初始轮廓线；④利用②所得到的梯度幅值定义snake模型的外力，实时计算snake曲线的能量E_m；将当前的snake曲线的能量E_m与上一时刻的snake曲线的能量E_m进行比较，当满足迭代参数时，所得到的snake曲线即为感兴趣物体的轮廓线。本发明依据得到的彩色图像梯度定义snake分割模型的图像外力，与传统的将彩色图像转为灰度图像处理或三通道分别处理的结果相比，能更快收敛到感兴趣物体的边缘，分割更准确。

Description

一种基于四元数的彩色snake图像分割方法

技术领域

本发明属于图像处理领域，涉及一种彩色图像分割方法，具体涉及一种基于四元数的彩色图像分割方法。

背景技术

随着近现代计算机和存储技术的飞速发展，各种电子数码照相设备得到普及，我们的日常生活都与各种电子数码影像技术相关，尤其是自然图像的处理得到越来越多的关注，而自然图像一般为彩色图像，彩色图像分割对于图像分析、理解和目标识别是至关重要的。

在传统的计算机视觉领域，严格各自独立的分层理论有着深远的影响，传统的理论认为底层的视觉任务只能依赖于由图像本身获得的信息来完成。1988年Kass等人首次提出主动轮廓模型(又被称作snake模型)，他们认为在许多的图像理解任务中，底层事件的正确理解依赖于高层知识。被称为snake的主动轮廓模型是能量极小化的样条，内力约束其形状，外力控制其变形行为，图像力将其引导至显著性的图像特征，snake模型的曲线能量将在图像特征附近极小化。Snake模型的成功之处还在于它对一系列广泛的视觉问题给出了统一的解决办法，从该模型首次提出，在近20年的研究中，snake模型已被成功应用于计算机视觉的许多领域，如边缘提取、图像分割和运动跟踪等。

然而该模型一直局限于灰度图像的处理，很少见诸于彩色图像的snake模型，究其主要原因是一直没有好的方法用来计算彩色图像的梯度，进而定义snake模型的图像外力。

发明内容

本发明的目的在于提供一种基于四元数的彩色snake图像分割方法，该方法将一个彩色像素作为一个整体进行处理，并在处理过程中很好地体现彩色三通道间的光谱联系。

本发明提供的一种基于四元数的彩色snake图像分割方法，其特征在于，该方法包括下述步骤：

第1步输入一幅待分割的RGB空间彩色图像，并将该图像表示到四元数空间；

第2步计算第1步中得到的四元数空间下表示的彩色图像的梯度幅值；

第3步在输入的待分割的RGB空间彩色图像上勾勒出感兴趣物体的初始轮廓线；

第4步利用第2步所得到的梯度幅值定义snake模型的外力，实时计算snake曲线的能量E_m；将当前的snake曲线的能量E_m与上一时刻的snake曲线的能量E_m进行比较，当满足迭代参数时，所得到的snake曲线即为感兴趣物体的轮廓线；

第2步中，按照下述过程计算得到梯度幅值：

用(R(x，y)，G(x，y)，B(x，y))表示待分割的RGB空间彩色图像，R(x，y)，G(x，y)，B(x，y)分别代表彩色图像的红绿蓝三通道，(x，y)代表像素坐标，则包含色度和亮度的彩色距离平方函数f(θ)写为：

f(θ)＝Ecos²θ+2Fcosθsinθ+Hsin²θ

其中θ(θ∈[-π/2，π/2])为图像梯度向量的角度方向，E，F，H定义如下：

$\left\{\begin{array}{c}E=E_1+E_2\\ F=F_1+F_2\\ H=H_1+H_2\end{array}\right.$

$\left\{\begin{array}{c}{E}_1=r_x^2+g_x^2+b_x^2\\ F_1=r_x{r}_y+g_x{g}_y+b_x{b}_y\\ H_1=r_y^2+g_y^2+b_y^2\end{array}\right.$

$\left\{\begin{array}{c}{E}_2={(w_1\frac{\∂R(x,y)}{\∂x}+w_2\frac{\∂G(x,y)}{\∂x}+w_3\frac{\∂B(x,y)}{\∂x})}^2\\ F_2=(w_1\frac{\∂R(x,y)}{\∂x}+w_2\frac{\∂G(x,y)}{\∂x}+w_3\frac{\∂B(x,y)}{\∂x})(w_1\frac{\∂R(x,y)}{\∂y}+w_2\frac{\∂G(x,y)}{\∂y}+w_3\frac{\∂B(x,y)}{\∂y})\\ H_2={(w_1\frac{\∂R(x,y)}{\∂y}+w_2\frac{\∂G(x,y)}{\∂y}+w_3\frac{\∂B(x,y)}{\∂y})}^2\end{array}\right.$

$\left\{\begin{array}{c}{r}_x=\frac{\∂R(x,y)}{\∂x}-\frac{1}{3}(\frac{\∂R(x,y)}{\∂x}+\frac{\∂G(x,y)}{\∂x}+\frac{\∂B(x,y)}{\∂x})\\ r_y=\frac{\∂R(x,y)}{\∂y}-\frac{1}{3}(\frac{\∂R(x,y)}{\∂y}+\frac{\∂G(x,y)}{\∂y}+\frac{\∂B(x,y)}{\∂y})\\ g_x=\frac{\∂G(x,y)}{\∂x}-\frac{1}{3}(\frac{\∂R(x,y)}{\∂x}+\frac{\∂G(x,y)}{\∂x}+\frac{\∂B(x,y)}{\∂x})\\ g_y=\frac{\∂G(x,y)}{\∂y}-\frac{1}{3}(\frac{\∂R(x,y)}{\∂y}+\frac{\∂G(x,y)}{\∂y}+\frac{\∂B(x,y)}{\∂y})\\ b_x=\frac{\∂B(x,y)}{\∂x}-\frac{1}{3}(\frac{\∂R(x,y)}{\∂x}+\frac{\∂G(x,y)}{\∂x}+\frac{\∂B(x,y)}{\∂x})\\ b_y=\frac{\∂B(x,y)}{\∂y}-\frac{1}{3}(\frac{\∂R(x,y)}{\∂y}+\frac{\∂G(x,y)}{\∂y}+\frac{\∂B(x,y)}{\∂y})\end{array}\right.$

其中E，F，H，E₁，F₁，H₁，E₂，F₂，H₂，r_x，g_x，b_x，r_y，g_y，b_y为临时变量，其中E，E₁，H₁，E₂，H₂∈[0，∞)；

则彩色距离平方函数f(θ)的最大值f_max和使得f(θ)取得最大值f_max的方向θ_max为：

$f_{\max }=\frac{1}{2}((E+H)+\sqrt{{(E-H)}^2+{\left(2F\right)}^2})$

${\mathrm{\θ}}_{\max }=\left\{\begin{array}{cc}\mathrm{sgn}\left(F\right)\arcsin {\left(\frac{f_{\max }-E}{{2f}_{\max }-E-H}\right)}^{1/2}+\mathrm{n\π},& {(E-H)}^2+F^2\≠0\\ \mathrm{Undefined},& {(E-H)}^2+F^2=0\end{array}\right.$

式中n取整数，其中sgn(□)为符号函数：

$\mathrm{sgn}\left(F\right)=\left\{\begin{array}{cc}1,& F\&GreaterEqual;0\\ -1,& F<0\end{array}\right.$

同时θ_max满足如下的约束关系：

${\mathrm{\&theta;}}_{\max }\&Element;\mathrm{n\&pi;}+\left\{\begin{array}{c}\left[\begin{array}{cc}0,& \mathrm{\&pi;}/4\end{array}\right],F\&GreaterEqual;0\mathrm{andE}-H\&GreaterEqual;0\\ \left[\begin{array}{cc}\mathrm{\&pi;}/4,& \mathrm{\&pi;}/2\end{array}\right],F\&GreaterEqual;0\mathrm{andE}-H<0\\ \left[\begin{array}{cc}-\mathrm{\&pi;}/2,& -\mathrm{\&pi;}/4\end{array}\right],F<0\mathrm{andE}-H<0\\ \left[\begin{array}{cc}-\mathrm{\&pi;}/4,& 0\end{array}\right],F<0\mathrm{andE}-H\&GreaterEqual;0\end{array}\right.$

f_max为所得的彩色图像梯度的幅值，θ_max为所得的彩色图像梯度的角度。

作为上述技术方案的改进，第3步具体包括下述过程：

设待分割的RGB空间彩色图像的大小为row×col×3，以该RGB空间彩色图像的几何中心为中心，以(row-1)和(col-1)作为矩形两邻边，在该RGB空间彩色图像上勾勒出一个包含感兴趣物体的矩形轮廓即为初始snake轮廓线v(s)＝(x(s)，y(s))，其中(x，y)代表构成初始轮廓线v(s)的像素的坐标，s∈(0，1)是将曲线进行参数化表示时的归一化参数。

作为上述技术方案的进一步改进，第4步具体包括下述过程：

第4.1步设定给定snake模型的形变参数，包括弹性系数α、刚性系数β和阻尼系数γ的值，α，β，γ∈[0，∞)，将v(s)视作s和时间t的函数，即将snake曲线v(s)视作一条动态曲线v(s，t)，给定时间步长1，迭代次数上限M，M为整数，且M∈[100，1500]，记初始勾勒snake曲线的时刻t＝t₀＝0，此时snake曲线能量记为E_snake＝E₀，令时间步的序号为m，其初始值为0；

第4.2步令m＝m+1，依据式一得到t_m时刻的snake曲线，再依据式二计算t_m时刻snake曲线的能量E_m：

$\left\{\begin{array}{c}{x}_{t_m}={(A+\mathrm{\&gamma;I})}^{-1}({\mathrm{\&gamma;x}}_{t_{m-1}}-F_x(x_{t_{m-1}},y_{t_{m-1}}))\\ y_{t_m}={(A+\mathrm{\&gamma;I})}^{-1}({\mathrm{\&gamma;y}}_{t_{m-1}}-F_y(x_{t_{m-1}},y_{t_{m-1}}))\end{array}\right.$ 式一

$E_m={\&Integral;}_0^1(E_{\mathrm{int}}\left(v\left(s\right)\right)+E_{\mathrm{ext}}\left(v\left(s\right)\right))\mathrm{ds}$ 式二

其中，E_int(v(s))表示snake曲线的内部能量，E_ext(v(s))表示snake曲线的外部能量，令E_ext(v(s))＝-f_max；

第4.3步将t_m时刻的snake曲线的能量E_m与t_m-1时刻的snake曲线能量E_m-1进行比较，若E_m＜E_m-1，则用t_m时刻的snake曲线更新替换t_m-1时刻的snake曲线，再进入第4.4步，否则直接执行第4.4步；

第4.4步，判断m≤M是否成立，如果是，则返回第4.2步，否则输出第M次循环所得snake曲线即为算法所得感兴趣物体的轮廓线。

本发明提供一种基于四元数的彩色snake图像分割方法。该方法首先利用四元数来表示彩色图像，可以将彩色像素当作一个整体进行处理，能很好地体现彩色三通道间的光谱联系。再依据本发明提出的一种基于四元数的彩色图像梯度计算方法来计算彩色图像的梯度，并以该梯度来定义snake模型的图像力，引导snake形变收敛到感兴趣物体的边缘，以达到分割出感兴趣物体的目的。

附图说明

图1为本发明一种基于四元数的彩色snake图像分割方法流程图；

图2为本发明一种基于四元数的彩色snake图像分割方法应用于一幅常见的彩色图像分割的实例；

图2a为常见的图像分割数据库中的一幅彩色图像；

图2b为应用本发明提出的基于四元数的彩色图像梯度计算方法得到的图2a的梯度幅值图；

图2c中沿着图像边缘的矩形框是本发明中由计算机自动确定的初始snake轮廓线；

图2d为本发明迭代一种基于四元数的彩色snake图像分割方法迭代1000次后得到的轮廓线。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明。在此需要说明的是，对于这些实施方式的说明用于帮助理解本发明，但并不构成对本发明的限定。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示，本发明提供的一种基于四元数的彩色snake图像分割方法，包括以下步骤：

第1步输入一幅待分割的RGB空间彩色图像R(x，y)，G(x，y)，B(x，y)，记图像大小为row×col×3，并将该图像表示到四元数空间。

四元数最初由Hamilton于1843年发现，四元数其实是复数的扩展。一个四元数(q)包含一个实部(q₀)和三个虚部(q₁，q₂，q₃)，经常被表示为如下形式：

q＝q₀+q₁i+q₂j+q₃k            (1)

其中q₀，q₁，q₂和q₃都是实数，i，j和k是三个虚部单位，并满足如下的运算法则：

$\left\{\begin{array}{c}{i}^2=j^2=k^2=\mathrm{ijk}=-1\\ \mathrm{ij}=k,\mathrm{jk}=i,\mathrm{ki}=j\\ \mathrm{ji}=-k,\mathrm{kj}=-i,\mathrm{ik}=-j\end{array}\right.---\left(2\right)$

四元数的模和共轭定义如下：

$\left|q\right|=\sqrt{q_0^2+q_1^2+q_2^2+q_3^2}---\left(3\right)$

q^*＝q₀-q₁i-q₂j-q₃k           (4)

如果一个四元数的实部为0，我们称该四元数为纯虚四元数；如果一个四元数的模为1，我们称该四元数为单位四元数；如果一个纯虚四元数的模为1，我们称该四元数为单位纯虚四元数。

RGB空间的一个彩色像素[r，g，b]^T在四元数空间可以用一个纯虚四元数表示：

q＝ri+gj+bk                  (5)

其中r，g，b分别对应RGB彩色空间中一个彩色像素的红绿蓝三通道的颜色值。

第2步计算第1步中得到的四元数空间下表示的彩色图像的梯度。

本发明提出一种全新的基于四元数的彩色图像梯度计算方法。

我们定义一个单位纯虚四元数

依据式(4)的定义，其共轭四元数为

依据式(2)的四元数运算法则，则UqU^*代表着在RGB三维颜色空间将一个彩色像素q围绕着灰度轴μ＝i+j+k(因为μ所代表的彩色像素红绿蓝三通道的颜色值相同，都为1，所以可以用颜色向量μ＝i+j+k代表灰度轴)旋转180度，这也就是说像素q和UqU^*在彩色空间中是关于灰度轴μ对称的。因此，依据向量加法的平行四边形法则和颜色相消原理，像素q+UqU^*应该落在灰度轴上。

我们现在假设有两个彩色像素分别用两个纯虚四元数表示为q₁＝r₁i+g₁j+b₁k和q₂＝r₂i+g₂j+b₂k，我们定义：

q₃＝q₂+Uq₁U^*＝r₃i+g₃j+b₃k                   (6)

那么我们可以利用q₃来定义如下的四元数用以表示两个像素q₁和q₂之间的色调差别：

$Q(q_1,q_2)=(r_3-\frac{r_3+g_3+b_3}{3})i+(g_3-\frac{r_3+g_3+b_3}{3})j+(b_3-\frac{r_3+g_3+b_3}{3})k---\left(7\right)$

通过简单的计算转换，上式可写为：

$Q(q_1,q_2)=((r_2-r_1)-\frac{1}{3}((r_2-r_1)+(g_2-g_1)+(b_2-b_1)))i+$

$((g_2-g_1)-\frac{1}{3}((r_2-r_1)+(g_2-g_1)+(b_2-b_1)))j+---\left(8\right)$

$((b_2-b_1)-\frac{1}{3}((r_2-r_1)+(g_2-g_1)+(b_2-b_1)))k$

我们再作如下定义：

CD(q₁，q₂)＝|Q(q₁，q₂)|²+|I(q₁，q₂)|²        (9)

I(q₁，q₂)＝w₁(r₂-r₁)+w₂(g₂-g₁)+w₃(b₂-b₁)        (10)

式(9)中CD(q₁，q₂)表示两个彩色像素q₁，q₂之间的颜色距离，该颜色距离由两部分组成：色度距离|Q(q₁，q₂)|²和亮度距离|I(q₁，q₂)|²，其中色度距离|Q(q₁，q₂)|²为式(8)所得四元数的模的平方，亮度距离为式(10)的绝对值的平方。式(10)中w_i(i＝1，2，3)，w_i∈(0，1)分别代表着红、绿、蓝三通道对图像亮度的贡献权值，典型情况下，我们取w₁＝w₂＝w₃＝1/3，意味着彩色图像亮度值为三个通道的平均值；在YIQ、YUV和YCbCr颜色空间中其亮度通道的值是由红绿蓝三通道值按照w₁＝0.299，w₂＝0.587，w₃＝0.114的比例计算所得。

在上面分析的基础上进一步研究，本发明提出一个全新的基于四元数的彩色图像梯度计算方法：

我们用(R(x，y)，G(x，y)，B(x，y))代表RGB彩色空间中的一幅彩色图像，R(x，y)，G(x，y)，B(x，y)分别代表彩色图像的红绿蓝三通道，(x，y)代表像素坐标，则包含色度和亮度的彩色距离平方函数f(θ)可写为：

f(θ)＝Ecos²θ+2Fcosθsinθ+Hsin²θ             (11)

其中θ(θ∈[-π/2，π/2])为图像梯度向量的角度方向，式(11)中E，F，H定义如下：

$\left\{\begin{array}{c}E=E_1+E_2\\ F=F_1+F_2\\ H=H_1+H_2\end{array}\right.---\left(12\right)$

$\left\{\begin{array}{c}{E}_1=r_x^2+g_x^2+b_x^2\\ F_1=r_x{r}_y+g_x{g}_y+b_x{b}_y\\ H_1=r_y^2+g_y^2+b_y^2\end{array}\right.---\left(13\right)$

$\left\{\begin{array}{c}{E}_2={(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;x}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;x}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;x})}^2\\ F_2=(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;x}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;x}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;x})(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;y}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;y}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;y})\\ H_2={(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;y}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;y}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;y})}^2\end{array}\right.---\left(14\right)$

$\left\{\begin{array}{c}{r}_x=\frac{\&PartialD;R(x,y)}{\&PartialD;x}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;x}+\frac{\&PartialD;G(x,y)}{\&PartialD;x}+\frac{\&PartialD;B(x,y)}{\&PartialD;x})\\ r_y=\frac{\&PartialD;R(x,y)}{\&PartialD;y}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;y}+\frac{\&PartialD;G(x,y)}{\&PartialD;y}+\frac{\&PartialD;B(x,y)}{\&PartialD;y})\\ g_x=\frac{\&PartialD;G(x,y)}{\&PartialD;x}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;x}+\frac{\&PartialD;G(x,y)}{\&PartialD;x}+\frac{\&PartialD;B(x,y)}{\&PartialD;x})\\ g_y=\frac{\&PartialD;G(x,y)}{\&PartialD;y}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;y}+\frac{\&PartialD;G(x,y)}{\&PartialD;y}+\frac{\&PartialD;B(x,y)}{\&PartialD;y})\\ b_x=\frac{\&PartialD;B(x,y)}{\&PartialD;x}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;x}+\frac{\&PartialD;G(x,y)}{\&PartialD;x}+\frac{\&PartialD;B(x,y)}{\&PartialD;x})\\ b_y=\frac{\&PartialD;B(x,y)}{\&PartialD;y}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;y}+\frac{\&PartialD;G(x,y)}{\&PartialD;y}+\frac{\&PartialD;B(x,y)}{\&PartialD;y})\end{array}\right.---\left(15\right)$

式(12)到(15)中E，F，H，E₁，F₁，H₁，E₂，F₂，H₂，r_x，g_x，b_x，r_y，g_y，b_y为临时变量，其中E，E₁，H₁，E₂，H₂∈[0，∞)。

则彩色距离平方函数f(θ)的最大值f_max和使得f(θ)取得最大值f_max的方向θ_max为：

$f_{\max }=\frac{1}{2}((E+H)+\sqrt{{(E-H)}^2+{\left(2F\right)}^2})---\left(16\right)$

${\mathrm{\&theta;}}_{\max }=\left\{\begin{array}{cc}\mathrm{sgn}\left(F\right)\arcsin {\left(\frac{f_{\max }-E}{{2f}_{\max }-E-H}\right)}^{1/2}+\mathrm{n\&pi;},& {(E-H)}^2+F^2\&NotEqual;0\\ \mathrm{Undefined},& {(E-H)}^2+F^2=0\end{array}\right.---\left(17\right)$

$\mathrm{sgn}\left(F\right)=\left\{\begin{array}{cc}1,& F\&GreaterEqual;0\\ -1,& F<0\end{array}\right.---\left(18\right)$

同时θ_max满足如下的约束关系：

${\mathrm{\&theta;}}_{\max }\&Element;\mathrm{n\&pi;}+\left\{\begin{array}{c}\left[\begin{array}{cc}0,& \mathrm{\&pi;}/4\end{array}\right],F\&GreaterEqual;0\mathrm{andE}-H\&GreaterEqual;0\\ \left[\begin{array}{cc}\mathrm{\&pi;}/4,& \mathrm{\&pi;}/2\end{array}\right],F\&GreaterEqual;0\mathrm{andE}-H<0\\ \left[\begin{array}{cc}-\mathrm{\&pi;}/2,& -\mathrm{\&pi;}/4\end{array}\right],F<0\mathrm{andE}-H<0\\ \left[\begin{array}{cc}-\mathrm{\&pi;}/4,& 0\end{array}\right],F<0\mathrm{andE}-H\&GreaterEqual;0\end{array}\right.---\left(19\right)$

式(17)和式(19)中，n取任一整数，其中sgn(□)为符号函数。

f_max即为本发明提出基于四元数的彩色梯度计算方法所得的彩色图像梯度的幅值，θ_max即为本发明提出基于四元数的彩色图像梯度计算方法所得的彩色图像梯度的角度。

第3步第1步中输入的彩色图像的大小为row×col×3，以输入图像的几何中心为中心，以(row-1)和(col-1)作为矩形两邻边，在输入图像上勾勒出一个包含感兴趣物体的矩形轮廓即为初始snake轮廓线v(s)＝(x(s)，y(s))，其中(x，y)代表构成初始轮廓线v(s)的像素的坐标，s∈(0，1)是将曲线进行参数化表示时的归一化参数。给定snake模型的形变参数(包括弹性系数α、刚性系数β和阻尼系数γ)的值，α，β，γ∈[0，∞)，将v(s)视作s和时间t的函数，即将snake曲线v(s)视作一条动态曲线v(s，t)，给定时间步长1，迭代次数上限M(M∈[100，1500]，M∈Z)，记初始勾勒snake曲线的时刻t＝t₀＝0，此时snake曲线能量记为E_snake＝E0，令时间步的序号为m，其初始值为0。

第4步令m＝m+1，依据式(32)得到t_m时刻的snake曲线，再依据式(33)计算t_m时刻snake曲线的能量E_m。

为了使式(33)中定义的曲线能量E_snake最小，由变分原理可知，曲线v(s)必须满足如下欧拉方程：

${\mathrm{\&alpha;v}}^{\&prime;\&prime;}\left(s\right)-\mathrm{\&beta;}{v}^{\&prime;\&prime;\&prime;\&prime;}\left(s\right)-\&dtri;E_{\mathrm{ext}}=0---\left(20\right)$

从物理学角度对snake模型进行分析，可以把式(20)看成为如下的力平衡方程：

F_int(v)+F_ext(v)＝0                            (21)

其中内力项F_int(v)定义为：

F_int(v)＝αv″(s)-βv″″(s)                  (22)

外力项F_ext(v)定义为：

$F_{\mathrm{ext}}\left(v\right)=-\&dtri;E_{\mathrm{ext}}---\left(23\right)$

为求解式(20)，我们将v(s)视作s和时间t的函数，即把曲线v(s)看作一个动态的曲线v(s，t)，根据牛顿第二定律，使轮廓线v(s，t)动态变化必须满足：

$m\frac{{\&PartialD;}^{2}v}{\&PartialD;t^2}=F_{\mathrm{int}}\left(v\right)+F_{\mathrm{ext}}\left(v\right)+F_{\mathrm{damp}}\left(v\right)---\left(24\right)$

其中F_damp(v)为阻尼力：

$F_{\mathrm{damp}}\left(v\right)=-\mathrm{\&gamma;}\frac{\&PartialD;v}{\&PartialD;t}---\left(25\right)$

m为质量单位，γ为阻尼系数。在图像分割中，为了防止惯性项

驱动轮廓超越边界，一般将质量单位m取为0，这样将式(25)代入式(24)得到snake的动力学方程为：

$\mathrm{\&gamma;}\frac{\&PartialD;v}{\&PartialD;t}=F_{\mathrm{int}}\left(v\right)+F_{\mathrm{ext}}\left(v\right)---\left(26\right)$

该方程即为snake曲线随时间变化的力平衡方程，当曲线的内力和外力相等时，曲线处于平衡状态停止运动，也即snake曲线达到能量最小化状态。

利用有限差分对式(20)进行离散化，可得到相应的欧拉方程：

α(v_i-v_i-1)-α(v_i+1-v_i)+β(v_i-2-2v_i-1+v_i)-2β(v_i-1-2v_i+v_i+1)

+β(v_i-2v_i+1+v_i+2)+(F_x(i)+F_y(i))＝0            (27)

其中：

$F_x\left(i\right)=\frac{{\&PartialD;E}_{\mathrm{ext}}}{\&PartialD;x}{|}_{x=x_i},$ $F_y\left(i\right)=\frac{{\&PartialD;E}_{\mathrm{ext}}}{\&PartialD;y}{|}_{y=y_i},$ i＝0，1...N-1                                      (28)

(注：snake曲线v(s)由N个像素点组成分别表示为v_i，v_i＝v(x_i(s)，y_i(s))，表示对梯度幅值图求x方向偏导数在v_i点取得的值，即x＝x_i的时候梯度幅值图x方向偏导数的值)。

N为组成snake曲线的像素点个数。我们可以将式(27)以矩阵的形式表示为如下形式：

$\left\{\begin{array}{c}\mathrm{Ax}+F_x=0\\ \mathrm{Ay}+F_y=0\end{array}\right.---\left(29\right)$

其中A_N×N为系数α和β所确定的稀疏带状对角矩阵，记

$\left\{\begin{array}{c}a=2\mathrm{\&alpha;}+6\mathrm{\&beta;}\\ b=-(\mathrm{\&alpha;}+4\mathrm{\&beta;})\\ c=\mathrm{\&beta;}\end{array}\right.$

则A_N×N可写为：

根据无惯性动力学方程式(26)可得：

$\left\{\begin{array}{c}{\mathrm{Ax}}_{t_m}+F_x=-\mathrm{\&gamma;}(x_{t_{m-1}}-x_{t_{m-1}})\\ {\mathrm{Ay}}_{t_m}+F_y=-\mathrm{\&gamma;}(y_{t_{m-1}}-y_{t_{m-1}})\end{array}\right.---\left(31\right)$

对式(31)做进一步的求解可得最终迭代解为：

$\left\{\begin{array}{c}{x}_{t_m}={(A+\mathrm{\&gamma;I})}^{-1}({\mathrm{\&gamma;x}}_{t_{m-1}}-F_x(x_{t_{m-1}},y_{t_{m-1}}))\\ y_{t_m}={(A+\mathrm{\&gamma;I})}^{-1}({\mathrm{\&gamma;y}}_{t_{m-1}}-F_y(x_{t_{m-1}},y_{t_{m-1}}))\end{array}\right.---\left(32\right)$

其中γ为阻尼系数，t为时间变量，具体求解是为离散的迭代次数，F_x，F_y定义见式(28)，A是由α和β系数所确定的稀疏带状对角矩阵，定义见式(30)，I为实单位矩阵。

依据式(32)得到更新的snake曲线后，利用第2步中式(16)所得彩色图像梯度幅值来定义snake模型的外部能量来计算曲线的能量，snake曲线的能量函数为：

$E_m=E_{\mathrm{snake}}\left(v\left(s\right)\right)={\&Integral;}_0^1{E}_{\mathrm{snake}}\left(v\left(s\right)\right)\mathrm{ds}$

$={\&Integral;}_0^1(E_{\mathrm{int}}\left(v\left(s\right)\right)+E_{\mathrm{ext}}\left(v\left(s\right)\right))\mathrm{ds}---\left(33\right)$

其中E_int(v(s))表示曲线的内部能量，E_ext(v(s))表示曲线的外部能量。内部能量为snake曲线提供张力和刚力，来控制曲线的连续性和光滑性，其中内部能量E_int(v(s))定义如下：

$E_{\mathrm{int}}\left(v\left(s\right)\right)=\frac{1}{2}[\mathrm{\&alpha;}{\left|v^{\&prime;}\left(s\right)\right|}^2+\mathrm{\&beta;}{\left|v^{\&prime;\&prime;}\left(s\right)\right|}^2]---\left(34\right)$

其中：

$v^{\&prime;}\left(s\right)=\frac{\&PartialD;v\left(s\right)}{\&PartialD;s},$ $v^{\&prime;\&prime;}\left(s\right)=\frac{{\&PartialD;}^{2}v\left(s\right)}{\&PartialD;s^2}$

分别代表snake轮廓线对归一化参变量s的一阶导数和二阶导数。

式(34)中α为弹性系数，一阶导数项v′(s)表示弹性能量项；β为刚性系数，二阶导数项v″(s)表示刚性能量。弹性能量的作用是控制变形曲线的伸缩，保持曲线的连续性，类似弹力绳的特性，刚性能量的作用是控制曲线的弯曲，使曲线保持平滑。当权重系数α和β都不为0时，则轮廓曲线是一条光滑的曲线，若α为0曲线会出现间断；若β为0曲线会产生角点；若β取值较大，会使得snake曲线能量E_snake(v(s))取极小值时是一个闭合的圆或者一条直线段。

Snake曲线的外部能量函数是定义在整幅图像平面(R(x，y)，G(x，y)，B(x，y))上的标量函数，反应了图像的某些特性，如边缘等。对于一幅给定的图像(R(x，y)，G(x，y)，B(x，y))，snake的外部能量E_ext(v(s))将主动地把轮廓线吸引到感兴趣物体的边界上，外部能量E_ext(v(s))函数表达式如下：

E_ext(v(s))＝E_image(v(s))+E_con(v(s))            (35)

式(35)中E_image(v(s))称作由图像力产生的能量，它将驱动snake曲线向感兴趣物体的真实轮廓移动；E_con(v(s))是由人为施加的外部力产生的能量，其目的是使得曲线能更快更精确地收敛到真实边界，外部力一般作为可选项。当snake模型用于感兴趣物体的分割时，我们利用图像的梯度来定义图像力产生的能量，而不施加外部控制力，则

E_ext(v(s))＝E_image(v(s))＝-f_max                (36)

其中f_max为第2步中通过式(16)计算所得的彩色图像梯度幅值。

第5步将t_m时刻的snake曲线的能量E_m与t_m-1时刻的snake曲线能量E_m-1进行比较，若E_m＜E_m-1，则用t_m时刻的snake曲线更新替换t_m-1时刻的snake曲线，再执行第6步；否则直接执行第6步；

第6步，判断m≤M是否成立，如果是，则返回第4步，否则输出第M次循环所得snake曲线即为算法所得感兴趣物体的轮廓线。

实例：

下面结合附图2的实例具体说明本发明一种基于四元数的彩色snake图像分割方法的步骤：

第1步输入一幅待分割的彩色图像，如附图2a，图2a的图像大小为481*321*3，我们用一个大小为481*321的四元数矩阵将其表示到四元数空间。

第2步利用本发明提出的基于四元数的彩色图像梯度计算方法计算图2a中图像的梯度，得到的梯度幅值图如图2b。

第3步程序自动地以输入图像的几何中心为中心做一个大小为480*320的线框作为初始snake轮廓线，初始化参数取α＝0.65，β＝0.001，γ＝0.3并依据式(33)计算snake曲线的能量，迭代次数上限取M＝1000。令时间步的序号为m，其初始值为0。

第4步令m＝m+1，依据式(32)得到t_m时刻的snake曲线，再依据式(33)计算t_m时刻snake曲线的能量Em

第5步将t_m时刻的snake曲线的能量E_m与t_m-1时刻的snake曲线能量E_m-1进行比较，若E_m＜E_m-1，则用t_m时刻的snake曲线更新替换t_m-1时刻的snake曲线，再执行第6步；否则直接执行第6步；

第6步，判断m≤M(在此实例中我们取M＝1000)是否成立，如果是，则返回第4步，否则输出第M次循环所得snake曲线即为算法所得感兴趣物体的轮廓线，如图2d。

以上所述为本发明的较佳实施例而已，但本发明不应该局限于该实施例和附图所公开的内容。所以凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。

Claims (3)

1.一种基于四元数的彩色snake图像分割方法，其特征在于，该方法包括下述步骤：

第1步输入一幅待分割的RGB空间彩色图像，并将该图像表示到四元数空间；

第2步计算第1步中得到的四元数空间下表示的彩色图像的梯度幅值；

第3步在输入的待分割的RGB空间彩色图像上勾勒出感兴趣物体的初始轮廓线；

第4步利用第2步所得到的梯度幅值定义snake模型的外力，实时计算snake曲线的能量E_m；将当前的snake曲线的能量E_m与上一时刻的snake曲线的能量E_m进行比较，当满足迭代参数时，所得到的snake曲线即为感兴趣物体的轮廓线；

第2步中，按照下述过程计算得到梯度幅值：

用(R(x，y)，G(x，y)，B(x，y))表示待分割的RGB空间彩色图像，R(x，y)，G(x，y)，B(x，y)分别代表彩色图像的红绿蓝三通道，(x，y)代表像素坐标，则包含色度和亮度的彩色距离平方函数f(θ)写为：

f(θ)＝Ecos²θ+2Fcosθsinθ+Hsin²θ

其中θ(θ∈[-π/2，π/2])为图像梯度向量的角度方向，E，F，H定义如下：

$\left\{\begin{array}{c}E=E_1+E_2\\ F=F_1+F_2\\ H=H_1+H_2\end{array}\right.$

$\left\{\begin{array}{c}{E}_1=r_x^2+g_x^2+b_x^2\\ F_1=r_x{r}_y+g_x{g}_y+b_x{b}_y\\ H_1=r_y^2+g_y^2+b_y^2\end{array}\right.$

$\left\{\begin{array}{c}{E}_2={(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;x}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;x}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;x})}^2\\ F_2=(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;x}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;x}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;x})(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;y}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;y}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;y})\\ H_2={(w_1\frac{\&PartialD;R(x,y)}{\&PartialD;y}+w_2\frac{\&PartialD;G(x,y)}{\&PartialD;y}+w_3\frac{\&PartialD;B(x,y)}{\&PartialD;y})}^2\end{array}\right.$

$\left\{\begin{array}{c}{r}_x=\frac{\&PartialD;R(x,y)}{\&PartialD;x}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;x}+\frac{\&PartialD;G(x,y)}{\&PartialD;x}+\frac{\&PartialD;B(x,y)}{\&PartialD;x})\\ r_y=\frac{\&PartialD;R(x,y)}{\&PartialD;y}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;y}+\frac{\&PartialD;G(x,y)}{\&PartialD;y}+\frac{\&PartialD;B(x,y)}{\&PartialD;y})\\ g_x=\frac{\&PartialD;G(x,y)}{\&PartialD;x}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;x}+\frac{\&PartialD;G(x,y)}{\&PartialD;x}+\frac{\&PartialD;B(x,y)}{\&PartialD;x})\\ g_y=\frac{\&PartialD;G(x,y)}{\&PartialD;y}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;y}+\frac{\&PartialD;G(x,y)}{\&PartialD;y}+\frac{\&PartialD;B(x,y)}{\&PartialD;y})\\ b_x=\frac{\&PartialD;B(x,y)}{\&PartialD;x}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;x}+\frac{\&PartialD;G(x,y)}{\&PartialD;x}+\frac{\&PartialD;B(x,y)}{\&PartialD;x})\\ b_y=\frac{\&PartialD;B(x,y)}{\&PartialD;y}-\frac{1}{3}(\frac{\&PartialD;R(x,y)}{\&PartialD;y}+\frac{\&PartialD;G(x,y)}{\&PartialD;y}+\frac{\&PartialD;B(x,y)}{\&PartialD;y})\end{array}\right.$

其中E，F，H，E₁，F₁，H₁，E₂，F₂，H₂，r_x，g_x，b_x，r_y，g_y，b_y为临时变量，其中E，E₁，H₁，E₂，H₂∈[0，∞)；

则彩色距离平方函数f(θ)的最大值f_max和使得f(θ)取得最大值f_max的方向θ_max为：

$f_{\max }=\frac{1}{2}((E+H)+\sqrt{{(E-H)}^2+{\left(2F\right)}^2})$

${\mathrm{\&theta;}}_{\max }=\left\{\begin{array}{cc}\mathrm{sgn}\left(F\right)\arcsin {\left(\frac{f_{\max }-E}{{2f}_{\max }-E-H}\right)}^{1/2}+\mathrm{n\&pi;},& {(E-H)}^2+F^2\&NotEqual;0\\ \mathrm{Undefined},& {(E-H)}^2+F^2=0\end{array}\right.$

式中n取整数，其中sgn(□)为符号函数：

$\mathrm{sgn}\left(F\right)=\left\{\begin{array}{cc}1,& F\&GreaterEqual;0\\ -1,& F<0\end{array}\right.$

同时θ_max满足如下的约束关系：

${\mathrm{\&theta;}}_{\max }\&Element;\mathrm{n\&pi;}+\left\{\begin{array}{c}\left[\begin{array}{cc}0,& \mathrm{\&pi;}/4\end{array}\right],F\&GreaterEqual;0\mathrm{andE}-H\&GreaterEqual;0\\ \left[\begin{array}{cc}\mathrm{\&pi;}/4,& \mathrm{\&pi;}/2\end{array}\right],F\&GreaterEqual;0\mathrm{andE}-H<0\\ \left[\begin{array}{cc}-\mathrm{\&pi;}/2,& -\mathrm{\&pi;}/4\end{array}\right],F<0\mathrm{andE}-H<0\\ \left[\begin{array}{cc}-\mathrm{\&pi;}/4,& 0\end{array}\right],F<0\mathrm{andE}-H\&GreaterEqual;0\end{array}\right.$

f_max为所得的彩色图像梯度的幅值，θ_max为所得的彩色图像梯度的角度。

2.根据权利要求1所述的基于四元数的彩色snake图像分割方法，其特征在于，第3步具体包括下述过程：

设待分割的RGB空间彩色图像的大小为row×col×3，以该RGB空间彩色图像的几何中心为中心，以(row-1)和(col-1)作为矩形两邻边，在该RGB空间彩色图像上勾勒出一个包含感兴趣物体的矩形轮廓即为初始snake轮廓线v(s)＝(x(s)，y(s))，其中(x，y)代表构成初始轮廓线v(s)的像素的坐标，s∈(0，1)是将曲线进行参数化表示时的归一化参数。

3.根据权利要求2所述的基于四元数的彩色snake图像分割方法，其特征在于，第4步具体包括下述过程：

第4.1步设定snake模型的形变参数，包括弹性系数α、刚性系数β和阻尼系数γ的值，α，β，γ∈[0，∞)，将v(s)视作s和时间t的函数，即将snake曲线v(s)视作一条动态曲线v(s，t)，给定时间步长1，迭代次数上限M，M为整数，且M∈[100，1500]，记初始勾勒snake曲线的时刻t＝t₀＝0，此时snake曲线能量记为E_snake＝E₀，令时间步的序号为m，其初始值为0；

第4.2步令m＝m+1，依据式一得到t_m时刻的snake曲线，再依据式二计算t_m时刻snake曲线的能量E_m：

$\left\{\begin{array}{c}{x}_{t_m}={(A+\mathrm{\&gamma;I})}^{-1}({\mathrm{\&gamma;x}}_{t_{m-1}}-F_x(x_{t_{m-1}},y_{t_{m-1}}))\\ y_{t_m}={(A+\mathrm{\&gamma;I})}^{-1}({\mathrm{\&gamma;y}}_{t_{m-1}}-F_y(x_{t_{m-1}},y_{t_{m-1}}))\end{array}\right.$ 式一

$E_m={\&Integral;}_0^1(E_{\mathrm{int}}\left(v\left(s\right)\right)+E_{\mathrm{ext}}\left(v\left(s\right)\right))\mathrm{ds}$ 式二

其中，E_int(v(s))表示snake曲线的内部能量，E_ext(v(s))表示snake曲线的外部能量，令E_ext(v(s))＝-f_max；

第4.3步将t_m时刻的snake曲线的能量E_m与t_m-1时刻的snake曲线能量E_m-1进行比较，若E_m＜E_m-1，则用t_m时刻的snake曲线更新替换t_m-1时刻的snake曲线，再进入第4.4步，否则直接执行第4.4步；

第4.4步，判断m≤M是否成立，如果是，则返回第4.2步，否则输出第M次循环所得snake曲线即为算法所得感兴趣物体的轮廓线。

Images