CN1367455A - 基于样条思内克斯(Snakes)模型的中医舌象计算机分割方法 - Google Patents

Abstract

基于样条思内克斯(Snakes)模型的中医舌象计算机分割方法,是将中医舌图象中的舌体从背景中分割出来以便于随后的特征分析。本发明中采用数码相机采集舌图象,并将图象输入到计算机进行处理、传输等操作。方法特征包括:1.根据对舌体形状的统计分析,定义一个以矩形区域为外边界的舌体轮廓模板;2.采用灰度投影分析法,获得一个矩形区域确定了舌体的大致位置和大小;3.提出一种基于灰度投影与刚性模板的舌体轮廓初始化方法;4.在思内克斯(Snakes)模型的能量函数中加入了有关舌体轮廓的先验知识;5.采用Catmull－Rom样条思内克斯(Snakes)模型表示舌体轮廓,并采用已有的Greedy方法进行优化,在彩色舌图象中得到舌体轮廓,分割出舌体区域。该分割准确、实用。

Description

基于样条思内克斯(Snakes)模型的中医舌象的计算机分割方法

技术领域

本发明涉及计算机医学图像处理领域，设计了一种基于样条思内克斯(Snakes)模型的舌体分割方法，将中医舌图象中的舌体从背景中分割出来，以便于随后的特征分析。

背景技术

常用在图像处理领域与舌象分割相关技术如下：

舌象分割即舌图象舌体区域分割，是图象分割方法的一种具体应用。传统的图象分割采用边缘检测、区域聚类等算法，这些方法均是利用图象的低层特征，即图像数据的均匀性或突变性，而没有利用有关分割目标的先验知识，如目标的位置、大小、形状等。这些方法需要进行复杂的后处理才能确定舌体的完整轮廓。而且由于轮廓特征和图象质量的差异等原因，往往难以获得令人满意的分割结果。

思内克斯(Snakes)模型又称主动轮廓线模型(active contour model)，由Kass等人于1987年提出。思内克斯(Snakes)是一种将目标的轮廓模型与图像特征相匹配的方法，它利用物体的多边形表示，把目标公式化为一个适当的能量函数E，然后将E最小化得到所期望的轮廓线，从而分割出图像中的目标。思内克斯(Snakes)，是一种能量极小化的样条，在各种力和约束条件的作用下发生弹性形变，直到得到期望的轮廓。其能量一般由三种力生成：内力约束它的形状，外力引导它的行为，图像力将其拖向显著的图像特征；Snakes的轮廓线锁定在图像特征附近，准确地将其极小化。这样，在寻找显著的图像特征时，高层机制可能通过将图像特征推向一个适当的局部极值点而与模型进行交互。思内克斯(Snakes)的最大特点是结合了图象的低层特征与高层知识，用于图象分割时，在鲁棒性、精确度、实用性等方面均优于传统方法。

一般地，Snakes由几种不同的作用力同时控制，每一种力产生一个能量项。可表示为 $E_{\mathrm{snake}}\left(V\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{\mathrm{snake}}\left(v_i\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[E_{\mathrm{internal}}\left(v_i\right)+E_{\mathrm{image}}\left(v_i\right)+E_{\mathrm{external}}\left(v_i\right)]$

其中E_internal称作内部能量，不依赖于图象数据，只与轮廓的形状特征有关，用来计算轮廓形状的一些我们所关注的特性，例如连续性和平滑性；E_image称作图象能量，是几何模型与图象数据相互作用产生的，是针对图象中感兴趣的特征定义的，例如：边缘、线、区域、纹理等；E_external为外部能量，代表各种人为定义的约束条件。

样条Snakes将一组离散坐标点作为控制点，以此定义一个样条曲线，这些点给出了曲线的大致形状。样条即样条曲线，是用一系列多项式曲线段形成的连续、光顺的曲线。在数学上，样条是一组分段平滑函数，用来近似或内插数据点集。样条Snakes所表示的曲线是解析、可微的，而且能够描述局部细节。可以解析求解它的几何微分特性，便于能量函数的计算。

采用思内克斯(Snakes)模型进行图象分割的整个过程是一个迭代优化直到收敛的过程。在此过程中需要解决三个问题：

a)轮廓的初始化，b)能量函数的设计，c)优化策略。

轮廓的初始化要确定轮廓点的初始位置和目标的初始形状和大小，通常采用人机交互的方式手工绘制。

优化策略决定了迭代收敛的速度。优化就是通过计算初始估计轮廓的内部能量和外部能量，按照一定的规则不断地对该估计轮廓的位置、形状等属性进行更新，直到轮廓满足某些预定的约束条件为止，这是一个迭代的过程。在这里，可以针对轮廓上的每一个控制点定义一个一定大小的搜索窗，在每一个搜索窗内选择具有最小Snakes能量的点取代原来的控制点作为新的轮廓控制点，依此类推，直到满足迭代终止条件。

优化是Snakes求解的一个重要步骤。初始估计轮廓在优化过程中逐渐向待测轮廓逼近，优化结束时，估计轮廓收敛于期望的目标轮廓。常用的优化算法有：变分法、动态规划法、Greedy算法、模拟退火算法等，其中Greedy算法是应用最广泛的方法。

Greedy方法是一种改进的动态规划方法。它实现能量函数的局部最优解，时间复杂度仅有O(nm)。因此计算速度快，且存储需求较小。对每一个点v_i，其一次迭代步骤如下：

1)计算E_snake(v_i)；

2)以v_i点m×m邻域内的其它点p_jk，j，k∈[1，m]取代v_i点，分别计算E_snake(p_jk)，j，k∈[1，m]；

3)记E_min＝min(E_snake(p_jk)，j，k∈[1，m])，若E_snake(v_i)＜E_min，则v_i点不变；否则将v_i点移至对应于E_min的邻域点。

迭代的终止准则为：不再产生控制点的移动或迭代次数超过预设的最大值。图6为Greedy算法的示意图，图中v_i为原点，v_i′为对应于E_min的邻域点。图中采用9邻域。

发明内容

本发明就是基于上述技术的研究，基于灰度投影与刚性模板的舌体轮廓自动初始化方法，实现了舌图象的自动分割。

为了将舌体区域自动、准确、快速地分割出来，本发明设计了一种基于思内克斯(Snakes)模型的舌体区域分割方法。这种方法考虑了舌体的形状、位置、大小、以及舌体区域与背景区的颜色差异等特点，设计能量函数，进行轮廓初始化，并且采用简便、快速的局部优化算法以达到实用。本发明的技术思路特征在于：

1、根据对舌体形状的统计分析，定义一个以矩形区域为外边界的舌体轮廓模板。通过引入该舌体模板，在分割时加入了关于舌体形状的先验知识，可以使轮廓曲线初始化时尽可能靠近真实的舌体轮廓。

2、采用灰度投影分析法，获得一个矩形区域，该区域确定了舌体的大位置和大小。

3、提出一种基于灰度投影与刚性模板的舌体轮廓初始化方法。

4、在思内克斯(Snakes)模型的能量函数中加入了有关舌体轮廓的先验知识。

5、经过舌体轮廓初始化并定义了能量函数后，采用Catmull-Rom样条思内克斯(Snakes)模型表示舌体轮廓，并采用已有的Greedy方法进行优化，从而在彩色舌图象中得到舌体轮廓，分割出舌体区域。

本发明的技术方案参见图1、图2。这种基于样条思内克斯(Snakes)模型的中医舌象的计算机分割方法，是由数码相机完成采集舌图象，并将舌体及色标的光学信号转换为电信号图象输入到计算机进行处理、传输等操作，其特征在于计算机处理主要是通过USB接口软件、在样条思内克斯(Snakes)模型的基础上对舌图像进行读/写处理，处理后的舌图像或对舌图象进行分割处理后输出到缓存器，经显示器显示结果，它依次包括下述步骤：

1)计算机从USB接口读入舌图象信号，并保存在内存中。

2)对舌图象进行变换，以增强舌体与周围背景之间的对比度，并将彩色图象变为灰度图象，变换针对各像素进行，变换公式为： $I(x,y)=\frac{R(x,y)-G(x,y)}{|G(x,y)-B(x,y)|+1}$

式中R(x，y)、G(x，y)和B(x，y)为像素的原始红、绿、蓝三色值，I(x，y)为变换后的灰度值。

因为我们发现舌图象像素的三色色度值R、G、B有以下规律：皮肤上的G值大于舌头边缘的G值，舌头边缘上的G值通常和B值差不多或较大，而在皮肤上的G值则均大于B值；舌头和皮肤的R值均大于G值和B值。根据这些彩色分量值彼此间的相对特性，为增强舌体轮廓与背景的灰度差异，我们采用以上的舌图象的强度变换函数。经过这样的变换，舌体部分的灰度值较高，而舌体以外的部分灰度值相对较低。

3)进入基于灰度投影与刚性模板，对样条思内克斯(Snakes)模型轮廓控制点进行初始化的轮廓初始化子程序：即采用灰度投影分析法，获得一个矩形区域，该区域确定了舌体的大致位置和大小。先根据增强图象水平和垂直方向的灰度或亮度投影的特征，确定一矩形区域的边界，由此获得矩形区域的上下左右4个边界，从而确定了舌体区域的位置和大小，在确定了舌体区域矩形边界后，计算舌体轮廓模板中外边界的刚性形变参数λ，进而确定初始控制点和轮廓，从而完成样条思内克斯(Snakes)模型的自动初始化。

4)为突出舌体，方便后续迭代过程的图象能量计算，对舌图象阈值化并进行非线性变换，分为以下几步：

①对2)中的增强图象求灰度值的最大值maxI，平均值meanI，并计算最大差值dI＝maxI-meanI；

②以meanI为阈值，对增强图象阈值化，将小于meanI的灰度值置为0；

③对阈值化后的图象进行非线性变换，变换公式为： $U(x,y)={\left(\frac{I(x,y)-\mathrm{meanI}}{\mathrm{dI}}\right)}^{\mathrm{\γ}}$

式中的指数γ为非线性变换参数。

5)用Greedy迭代优化算法对样条思内克斯(Snakes)模型进行求解，求解步骤如前所述，直到满足终止条件：一次迭代后轮廓控制点的位置不再发生变化，或迭代次数达到某个预设的最大值，轮廓控制点迭代优化子程序调用结束，得到舌体轮廓控制点的最终位置，采用通用Catmull-Rom样条内插公式进行样条内插，即可得到连续的舌体轮廓曲线。

6)采用通用Catmull-Rom样条内插公式进行样条内插，得到连续轮廓；思内克斯(Snakes)的能量最小化是在稀疏的离散点上进行的，为了得到目标的连续轮廓，需要进行插值得到精确、连续的轮廓线。

7)按通常方法将位于舌体轮廓曲线上像素的三色色度值置为(255，255，255)(白色)，保存为结果文件并输出。

另根据以上所述的基于样条思内克斯(Snakes)模型的中医舌象的计算机分割方法，其中对样条思内克斯(Snakes)模型轮廓控制点进行初始化的轮廓初始化子程序的特征在于，分五步进行轮廓初始化：

1)在增强图象的基础上，进行水平方向的灰度投影，利用舌体区域灰度投影值明显增高的特点，从投影的中部开始，分别向左、向右搜索第一个投影值突然减小的位置，获得矩形区域的左右边界；

2)左右边界的中心定为水平中线，将中线附近的局部区域定义为中心区；

3)将中心区进行垂直方向的灰度投影，由于舌体区下边界附近的增强灰度值常发生明显变化，因此下边界对应于投影值的一阶导数最大值；

4)将中心区进行垂直方向的原始图象亮度值(即(R+G+B)/3)投影，由于存在唇下阴影，舌体区域的上边界对应于投影值的最小值；

5)以上步骤确定矩形边界后，根据舌体轮廓模板中外边界的刚性形变参数，确定轮廓控制点的初始位置，形变参数的计算方法为：设刚性模板的外边界矩形长为L₀，宽为W₀，由以上四步得到的矩形左右上下边界的在舌图象中的位置分别为：l₁、r₁、t₁、b₁，则矩形区域长为L＝r-l，宽为W＝b-t，那么形变参数为：λ＝(λ_l，λ_w)，λ₁＝L/L₀，λ_w＝W/W₀，由于轮廓模板的控制点位置是相对于外边界矩形定义，即以该矩形的左上角顶点为坐标中心(0，0)，所以若设模板中第i个轮廓控制点的坐标位置为(x_0i，y_0i)，则得到的初始化轮廓控制点的坐标为：

                    x_i＝l+λ_lx_0i

                    y_i＝t+λ_wy_0i。

另根据以上所述的基于样条思内克斯(Snakes)模型的中医舌象的计算机分割方法的特征在于：内部能量采用通用形式，在计算图象能量时，首先对原始图象阈值化并进行非线性变换以突出舌体，然后根据阈值化后的图象强度和梯度计算图象能量，内部能量与图象能量之和为思内克斯(Snakes)模型的总能量。

对于某一轮廓控制点v_i，考虑到舌图象中舌体轮廓符合连续性与平滑性的假设，内部能量采用通用形式，即：E_internal(v_i)＝α(v_i)E_elas(v_i)+β(v_i)E_bend(v_i) $E_{\mathrm{elas}}\left(v_i\right)=|\frac{\∂V}{\∂s}{|}^2\≈\frac{|v_i-v_{i-1}{|}^2}{{\mathrm{ds}}^2}=\frac{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}{{\mathrm{ds}}^2}$ $E_{\mathrm{bend}}\left(v_i\right)=|\frac{{\∂}^{2}V}{\∂s^2}{|}^2\≈\frac{|v_{i-1}-2v_i+v_{i+1}{|}^2}{{\mathrm{ds}}^4}=\frac{{(x_{i-1}-2x_i+x_{i+1})}^2+{(y_{i-1}-2y_i+y_{i+1})}^2}{{\mathrm{ds}}^4}$

其中权值参数α、β分别控制对连续性和平滑性的约束程度。

图象能量为E_image(v_i)＝γ₁U(v_i)+γ₂|U(v_i))|²

其中权值参数γ₁、γ₂控制对轮廓所在区域的I(x，y)强度约束和梯度约束；γ₁、γ₂设为负值，使能量最小值趋近于的极大值以及梯度响应的极大值，

外部能量E_external＝0，

总能量函数表示为 $E\left(v_i\right)=\mathrm{\α}\frac{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}{{\mathrm{ds}}^2}+\mathrm{\β}\frac{{(x_{i-1}-2x_i+x_{i+1})}^2+{(y_{i-1}-2y_i+y_{i+1})}^2}{{\mathrm{ds}}^4}$ +γ₁U(v_i)+γ₂|U(v_i))|

另说明2点：

1、图像处理中，为了突出具体图象的某些特征，常常要对图象做一些变换，即对图象像素的色度值或灰度值进行运算，又称这样的变换为图象增强。例如，舌图象像素的三色色度值R、G、B有以下规律：皮肤上的G值大于舌头边缘的G值，舌头边缘上的G值通常和B值差不多或较大，而在皮肤上的G值则均大于B值；舌头和皮肤的R值均大于G值和B值。根据这些彩色分量值彼此间的相对特性，为增强舌体轮廓与背景的灰度差异，我们采用舌图象的强度变换函数I(x，y)： $I(x,y)=\frac{R(x,y)-G(x,y)}{|G(x,y)-B(x,y)|+1}$

在计算能量函数时所进行的非线性变换也是基于同样的原因，引入指数变换是为了更加突出舌体和皮肤的对比度。

2、轮廓的初始化要确定轮廓点的初始位置和目标的初始形状和大小，通常采用人机交互的方式手工绘制。由于在本系统要分割舌图象中的舌体，如果够自动确定舌体的大致位置、大小和形状，则轮廓初始化可以自动进行。本发明就是基于这样的考虑，基于灰度投影与刚性模板的舌体轮廓自动初始化方法，实现了舌图象的自动分割。能量函数要根据图像和分割目标的特点进行设计。本发明的效果见图7、图8。可以看出采用基于样条思内克斯(Snakes)模型，并结合舌图象的特点，可以自动得到准确的舌体轮廓。不需手工初始化，从而为随后的精确的舌象分析创造了条件。

附图说明

图1是中医舌象分割系统框图。

1、数码相机，2、USB接口，3、计算机处理器，4、输出缓存，5、舌体分割，6、显示器，7、分割结果；

图2是本发明方法主程序流程图；

图3是本发明方法中轮廓初始化子程序流程图；

图4是本发明方法中迭代优化子程序流程图；

图6是本发明中Greedy迭代优化示意图；

图7是轮廓的初始化过程示例

图中(a)水平灰度投影、中线以及中心区示意；(b)水平灰度投影与左、右边界示意；(c)垂直灰度投影的一阶导数与下边界示意；(d)垂直彩色亮度投影与上边界示意；(e)初始化结果；

图8是本发明中一幅彩色舌图象及其分割结果(白线表示轮廓线)；

图9是本发明方法在计算机上运行的系统主程序流程图；

图10是本发明方法在计算机上运行的轮廓初始化子程序流程图；

图11是本发明方法在计算机上运行的迭代优化子程序流程图。

具体实施方式

在图1的中医舌象分割系统框图中，数码相机和USB接口都是市售的，主要完成采集舌图象，将舌体及色标的光学信号转换为图象电信号输入到计算机，便于计算机处理、传输等操作；计算机处理主要是通过USB接口软件对舌图像进行读/写处理；处理后的舌图像输出到缓存器，便于显示；显示器是图象的输出设备，人眼通过显示器观看原始图象和分割后图象；舌体分割是对计算机读入的舌图象进行分割处理，输出结果。原始舌图象可以是通过数码相机实时采集到的图象，也可以是实现通过数码相机采集后保存在计算机硬盘里的图象。舌体分割系统在计算机中完成以下主程序，参见图8、图9：1、在有色标的情况下：初始化参数l、r、t、b、δ、γ、α、β、γ₁、γ₂、w、N及舌体模板。其中l、r、t、b为预设的值，是为了去掉色标区9色标：图8周边的三个色块条)和上嘴唇以上的部分。它们的值可以事先根据一幅舌图象确定。δ为求左右边界所用的阈值。根据试探法事先设定。本发明中δ＝0.2。γ为非线性变换所采用的参数，此处设γ＝0.15。α、β、γ₁、γ₂为思内克斯(Snakes)的能量计算的加权系数，四个参数的取值由实验获得。在本发明中可分别设为：α＝0.5，β＝3.0，γ₁＝-0.2，γ₂＝-1.0。w为思内克斯(Snakes)迭代优化时窗口的宽度，本发明中设w＝3。N是迭代优化时设定的最多迭代次数，本发明中N＝200。舌体模板用一个表示点相对坐标的数组来表示。在无色标的情况下：初始化参数l＝0、t＝0，b为图象高，r为图象宽。

2、对舌图象进行变换。变换逐个像素进行。变换公式为： $I(x,y)=\frac{R(x,y)-G(x,y)}{|G(x,y)-B(x,y)|+1}.$

3、舌体轮廓初始化。分为以下各步：

第一步：计算增强图像的水平方向投影。如前所述，增强图象I(x，y)有h行w列。从第l列开始直到第r列，计算每一列上所有像素(共r-l+1个像素)的灰度值之和S_i，

得到一数组S，S＝[S₁，S₂...S_i...S_r]，该数组即表示增强图像的水平方向灰度投影。其中l、r、t、b为预设的值，是为了去掉色标区和上嘴唇以上的部分，简化计算并充分利用前述舌体的颜色特点，而不受色标区的影响。它们的值可以事先根据一幅舌图象确定。

第二步：利用舌体区域灰度投影明显增高的特点，获得矩形区域的左右边界。方法是：

从投影数组一半处j＝[n/2]开始，向左方搜索(即每次令j减小1)第一个小于δ_l的S_j。其中[n/2]表示对n/2取整。即如果

       S_j＜δ_i         (0＜j＜[n/2])

则直线x＝l₁，l₁＝i+l被定为矩形区域的左边界。

从投影数组一半处j＝[n/2]开始，向右方搜索(即每次令j增大1)第一个达于δ_r的S_j。即如果

      S_j＜δ_l    ([n/2]+1＜j＜r-l+1)

则直线x＝r₁，r₁＝j+l被定为矩形区域的左边界。

边界确定后，即停止在这一区域的搜索。式中，δ_l，δ_r为预设的阈值，在此设δ_l＝δ_r＝δ(maxS-minS)+minS。maxS、minS分别是水平投影的最大、最小值。δ事先设定。本发明中δ＝0.2。

第三步：将左右边界的中线x＝[(l₁+r₁)/2]定为水平中线，将中线附近的局部区域定义为中心区。设中心区域的宽度为w₁，即从[(l₁+r₁)/2]-w₁到[(l₁+r₁)/2]+w₁为中心区。中心区域的宽度可以定为左右边界差的一半。

第四步：确定矩形区域下边界。如前所述，将中心区进行垂直方向的灰度投影，由于舌体区下边界附近的增强灰度值常发生明显变化，因此下边界对应于投影值的一阶导数最大值。按上述第一步的方法计算中心区垂直方向投影。设垂直投影数组为P，P＝[P₁，P₂...P_i...P_r]，则一阶导数的计算方法为：

            P_i＝P_i-P_i-1，  0＜i＜h

求一阶导数的最大值，根据前述，一阶导数最大值所对应的位置即为下边界b₁。

第五步：将中心区进行垂直方向的原始图象亮度值(即(R+G+B)/3)投影V，V的最小值所对应的位置即为舌体区域的上边界t₁。

第六步：计算轮廓模板的形变参数，

            λ＝(λ₁，λ_w)，λ_l＝L/L₀，λ_w＝W/W₀

按下式计算轮廓控制点的初始位置：

                     x_i＝l+λ_lx_0i

                     y_i＝t+λ_wy_0i

(x_0i，y_0i)为模板中第i个轮廓控制点的坐标位置。在本发明中，共有32个轮廓控制点。该组数据保存在一个数组中。

图10为本发明方法在计算机上运行的轮廓初始化子程序流程图，图7为轮廓的初始化过程示例图。

4、对舌图象阈值化并进行非线性变换。分为以下几步：

1)对2中获得的增强图象求灰度值的最大值maxI，平均值meanI，并计算最大差值dI＝maxI-meanI

2)meanI为阈值，对增强图象阈值化，将小于meanI的灰度值置为0；

3)对阈值化后的图象进行非线性变换，变换公式为： $U(x,y)={\left(\frac{I(x,y)-\mathrm{meanI}}{\mathrm{dI}}\right)}^{\mathrm{\γ}}$

式中的指数γ为非线性变换参数，由试探法获得。在本发明中γ＝0.15。

5、对舌体轮廓迭代优化。

迭代过程见图11，描述如下：

1)读入轮廓控制点初始坐标。将迭代次数i置为0，将指示是否有控制点位置变动的变量bMove置为假；

2)读入一个轮廓控制点的坐标；

3)在大小为w×w的搜索窗内逐个计算控制点邻域内各点(共w×w个点)的能量，并求最小能量，并记录最小能量所对应的点的坐标。本发明中，搜索窗采用3邻域，即w＝3。在能量函数的计算公式： $E\left(v_i\right)=\mathrm{\α}\frac{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}{{\mathrm{ds}}^2}+\mathrm{\β}\frac{{(x_{i-1}-2x_i+x_{i+1})}^2+{(y_{i-1}-2y_i+y_{i+1})}^2}{{\mathrm{ds}}^4}$

+γ₁U(v_i)+γ₂|U(v_i))|中，ds＝1，四个参数的取值由实验获得，分别为：

        α＝0.5，β＝3.0，γ₁＝-0.2，γ₂＝-1.0 。

4)将该轮廓控制点移至能量最小点。即将该控制点的坐标修改为最小能量所对应的点的坐标。若修改后的坐标与原来的坐标(步骤二)中读入的坐标)不同，则令bMove为真，否则不改变bMove的值；

5)如果轮廓控制点已经处理完毕，则令迭代次数i增加1；否则转二)，处理下一个轮廓控制点；

6)根据bMove的真假判断是否有控制点的坐标发生变化。如果没有，则迭代结束，返回，否则；

7)判断迭代次数i是否小于预设的最大迭代次数N，如果i＜N，则转二)，进行下一轮迭代，否则迭代结束，返回。

6、根据控制点坐标进行样条插值获得连续的轮廓曲线，可以采用不同的样条。本发明中采用Catmull-Rom样条样条内插方法。

例如Catmull-Rom样条是一种局部内插样条，其第i段Catmull-Rom样条可表示为： $Q_i\left(\mathrm{\τ}\right)={\mathrm{\ω}}_i^1\left(\mathrm{\τ}\right)Q_i^1\left(\mathrm{\τ}\right)+{\mathrm{\ω}}_i^2\left(\mathrm{\τ}\right)Q_i^2\left(\mathrm{\τ}\right)$

其中ω_i ¹(τ)、ω_i ²(τ)为样条的线性弯曲函数，定义为： ${\mathrm{\ω}}_i^1\left(\mathrm{\τ}\right)={\{}_{0---------\mathrm{otherwise}}^{\frac{{\mathrm{\τ}}_{i+1}-\mathrm{\τ}}{{\mathrm{\τ}}_{i+1}-{\mathrm{\τ}}_i}---\mathrm{if}{\mathrm{\τ}}_i\≤\mathrm{\τ}\≤{\mathrm{\τ}}_{i+1}}$ ${\mathrm{\ω}}_i^2\left(\mathrm{\τ}\right)={\{}_{0---------\mathrm{otherwise}}^{\frac{\mathrm{\τ}-{\mathrm{\τ}}_i}{{\mathrm{\τ}}_{i+1}-{\mathrm{\τ}}_i}---\mathrm{if}{\mathrm{\τ}}_i\≤\mathrm{\τ}\≤{\mathrm{\τ}}_{i+1}}$

Q_i ¹(τ)为内插v_i-1，v_i，v_i+1的二阶多项式，Q_i ²(τ)为内插v_i-1，v_i，v_i+1的二阶多项式，式，分别定义为： $Q_i^1\left(\mathrm{\τ}\right)=v_{i-1}+(\mathrm{\τ}-{\mathrm{\τ}}_{i-1})(\frac{-v_{i-1}+v_i}{-{\mathrm{\τ}}_{i-1}+{\mathrm{\τ}}_i}+\frac{(\mathrm{\τ}-{\mathrm{\τ}}_i)(-\frac{-v_{i-1}+v_i}{-{\mathrm{\τ}}_{i-1}+{\mathrm{\τ}}_i}+\frac{-v_i+v_{i+1}}{-{\mathrm{\τ}}_i+{\mathrm{\τ}}_{i+1}})}{-{\mathrm{\τ}}_{i-1}+{\mathrm{\τ}}_{i+1}})$ $Q_i^2\left(\mathrm{\τ}\right)=v_i+(\mathrm{\τ}-{\mathrm{\τ}}_i)(\frac{-v_i+v_{i+1}}{-{\mathrm{\τ}}_i+{\mathrm{\τ}}_{i+1}}+\frac{(\mathrm{\τ}-{\mathrm{\τ}}_{i+1})(-\frac{-v_i+v_{i+1}}{-{\mathrm{\τ}}_i+{\mathrm{\τ}}_{i+1}}+\frac{-v_{i+1}+v_{i+2}}{-{\mathrm{\τ}}_{i+1}+{\mathrm{\τ}}_{i+2}})}{-{\mathrm{\τ}}_i+{\mathrm{\τ}}_{i+2}})$

样条模型的参数τ的定义采用非均匀参数化方法，使其正比于控制点间的Euclidean距离，递推公式为： $\frac{{\mathrm{\τ}}_i-{\mathrm{\τ}}_{i-1}}{{\mathrm{\τ}}_{i+1}{-\mathrm{\τ}}_i}\text{=||}\frac{v_i-v_{i-1}}{v_{i+1}{-v}_i}\left|\right|$ 可选τ₁＝0和τ_n＝1。并将最终得到的轮廓线用白线表示(即将该像素的三色值置为(255，255，255))。

7、输出分割结果

图8是一幅舌图象分割结果。

Claims (3)

1、基于样条思内克斯(Snakes)模型的中医舌象计算机分割方法，是由数码相机完成采集舌图象，并将舌体及色标的光学信号转换为图象电信号输入到计算机进行处理、传输等操作，该方法特征在于计算机处理主要是通过USB接口软件对舌图像进行读/写处理，对处理后的舌图像在样条思内克斯(Snakes)模型的基础上进行分割处理后输出到缓存器，经显示器显示结果，它依次包括下述步骤：

1)计算机从USB接口读入舌图象信号，并保存在内存中；

2)对舌图象进行变换，以增强舌体与周围背景之间的对比度，并将彩色图象变为灰度图象，变换针对各像素进行，变换公式为： $I(x,y)=\frac{R(x,y)-G(x,y)}{|G(x,y)-B(x,y)|+1}$

式中R(x，y)、G(x，y)和B(x，y)为像素的原始红、绿、蓝三色值，I(x，y)为变换后的灰度值；

3)进入基于灰度投影与刚性模板，对样条思内克斯(Snakes)模型轮廓控制点进行初始化的轮廓初始化子程序：即采用灰度投影分析法，获得一个矩形区域，该区域确定了舌体的大致位置和大小，先根据增强图象水平和垂直方向的灰度或亮度投影的特征，确定一矩形区域的边界，由此获得矩形区域的上下左右4个边界，从而确定了舌体区域的位置和大小，在确定了舌体区域矩形边界后，计算舌体轮廓模板中外边界的刚性形变参数λ，进而确定初始控制点和轮廓，从而完成样条思内克斯(Snakes)模型的自动初始化；

4)对舌图象阈值化并进行非线性变换，分为以下几步：

①对2)中的增强图象求灰度值的最大值maxI，平均值meanI，并计算最大差值dI＝maxI-meanI；

②以meanI为阈值，对增强图象阈值化，将小于meanI的灰度值置为0；

③对阈值化后的图象进行非线性变换，变换公式为： $U(x,y)={\left(\frac{I(x,y)-\mathrm{meanI}}{\mathrm{dI}}\right)}^{\mathrm{\γ}}$

式中的指数γ为非线性变换参数；

5)用Greedy迭代优化算法对样条思内克斯(Snakes)模型进行求解，求解步骤采用通用方法，直到满足终止条件，一次迭代后轮廓控制点的位置不再发生变化，或迭代次数达到某个预设的最大值，轮廓控制点迭代优化子程序调用结束，得到舌体轮廓控制点的最终位置，采用通用Catmull-Rom样条内插公式进行样条内插，即可得到连续的舌体轮廓曲线；

6)采用通用Catmull-Rom样条内插公式进行样条内插，得到连续轮廓；

7)按通常方法将位于舌体轮廓曲线上像素的三色色度值置为(255，255，255)(白色)，保存为结果文件并输出。

2、根据权利要求1所述的基于样条思内克斯(Snakes)模型的中医舌象计算机分割方法，其中对样条思内克斯(Snakes)模型轮廓控制点进行初始化的轮廓初始化子程序的特征在于，分五步进行轮廓初始化：

1)在增强图象的基础上，进行水平方向的灰度投影，利用舌体区域灰度投影明显增高的特点，从投影的中部开始，分别向左、向右搜索第一个投影突然减小的位置，获得矩形区域的左右边界；

2)左右边界的中心定为水平中线，将中线附近的局部区域定义为中心区；

3)将中心区进行垂直方向的灰度投影，由于舌体区下边界附近的增强灰度值常发生明显变化，因此下边界对应于投影值的一阶导数最大值；

4)将中心区进行垂直方向的原始图象亮度值(即(R+G+B)/3)投影，由于存在唇下阴影，舌体区域的上边界对应于投影值的最小值；

5)以上步骤确定矩形边界后，根据舌体轮廓模板中外边界的刚性形变参数，确定轮廓控制点的初始位置，形变参数的计算方法为：设刚性模板的外边界矩形长为L₀，宽为W₀，由以上四步得到的矩形左右上下边界的在舌图象中的位置分别为：l₁、r₁、t₁、b₁，则矩形区域长为L＝r-l，宽为W＝b-t，那么形变参数为：λ＝(λ_l，λ_w)，λ_l＝L/L₀，λ_w＝W/W₀，由于轮廓模板的控制点位置是相对于外边界矩形定义，即以该矩形的左上角顶点为坐标中心(0，0)，所以若设模板中第i个轮廓控制点的坐标位置为(x_0i，y_0i)，则得到的初始化轮廓控制点的坐标为：

                    x_i＝l+λ_lx_0i

                    y_i＝t+λ_wy_0i。

3、根据权利要求1所述的基于样条思内克斯(Snakes)模型的中医舌象计算机分割方法，其特征在于：内部能量采用通用形式，在计算图象能量时，首先对原始图象阈值化并进行非线性变换以突出舌体，然后根据阈值化后的图象强度和梯度计算图象能量，内部能量与图象能量之和为思内克斯(Snakes)模型的总能量；

对于某一轮廓控制点v_i，内部能量采用通用形式，即：E_internal(v_i)＝α(v_i)E_elal(v_i)+β(v_i)E_bend(v_i) $E_{\mathrm{elas}}\left(v_i\right)=|\frac{\∂V}{\∂s}{|}^2\≈\frac{|v_i-v_{i-1}{|}^2}{{\mathrm{ds}}^2}=\frac{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}{{\mathrm{ds}}^2}$ $E_{\mathrm{bend}}\left(v_i\right)=|\frac{{\∂}^{2}V}{\∂s^2}{|}^2\≈\frac{|v_{i-1}-2v_i+v_{i+1}{|}^2}{{\mathrm{ds}}^4}=\frac{{(x_{i-1}-2x_i+x_{i+1})}^2+{(y_{i-1}-2y_i+y_{i+1})}^2}{{\mathrm{ds}}^4}$

其中权值参数α、β分别控制对连续性和平滑性的约束程度，

图象能量为：

         E_image(v_i)＝γ₁U(v_i)+γ₂|U(v_i))|²

其中权值参数γ₁、γ₂控制对轮廓所在区域的I(x，y)强度约束和梯度约束，γ₁、γ₂设为负值，使能量最小值趋近于的极大值以及梯度响应的极大值，

外部能量E_external＝0

总能量函数表示为： $E\left(v_i\right)=\mathrm{\α}\frac{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}{{\mathrm{ds}}^2}+\mathrm{\β}\frac{{(x_{i-1}-2x_i+x_{i+1})}^2+{(y_{i-1}-2y_i+y_{i+1})}^2}{{\mathrm{ds}}^4}$ +γ₁U(v_i)+γ₂|U(v_i))|

Images