CN1920854A

CN1920854A - 中医舌诊用的舌头彩色数码照片的舌体分割提取方法

Info

Publication number: CN1920854A
Application number: CNA2006101131195A
Authority: CN
Inventors: 白净; 张永红; 吴佳
Original assignee: Tsinghua University
Current assignee: Tsinghua University
Priority date: 2006-09-15
Filing date: 2006-09-15
Publication date: 2007-02-28
Anticipated expiration: 2026-09-15
Also published as: CN100375979C

Abstract

本发明属于中医舌诊用的舌头彩色图像处理技术领域。具体特征在于，用两步对舌体进行分割提取：首先得到舌体的粗略轮廓图，其次用二维形变轮廓模型进行较准确的分割。在第一步中，用分水岭变换对图像的红色通道进行处理，得到粗略的舌体轮廓，在处理时先用饱和度图阈值求图像的标记图像，再用标记图像对图像的梯度图进行强制最小值，得到待分割图像，然后再用分水岭变换方法对待分割图像采用在图像纵剖面的采样点—灰度值曲线上，实施分水岭变换，得到一个由各纵剖面上的分水岭构成的集合，该分水岭集合中的一片连通区域即为舌体区域，其边界就是舌体粗略轮廓图。试验证明，图像结果非常准确，适用于不同年龄层次的健康和疾病人群。

Description

中医舌诊用的舌头彩色数码照片的舌体分割提取方法

技术领域

本发明属于中医舌诊用的舌头彩色图像处理技术领域。

背景技术

从20世纪80年代中后期开始，中医工作者和数字图像工作者开始联合，尝试以计算机技术为基础进行舌诊客观化研究。研究主要集中在舌头照片颜色的校准，存储和输出和用现代图像处理分析技术对舌头的照片进行处理等方面。其中，用图像处理技术对舌头的各个参数进行分析的前提条件就是从人脸舌头的图像中提取出舌体区域部分，把其他部分去除。针对这个课题，人们采用了很多方法进行了研究。

之前所采用的分割提取算法的不足之处主要为：

1)早期的分割算法采用简单的域值分割结合数学形态学的办法进行舌体的提取，分割的结果舌头边缘轮廓曲线比较粗糙，不能够精确的勾勒出舌头的边缘。

2)中后期采用二维形变轮廓模型进行舌体的分割，舌头轮廓边缘较准确。但该方法的初始曲线取得的办法有些为人体手工描绘，不能够全自动；有些为指定矩形区域，导致曲线迭代时间过长，方法太过于耗时。

3)鉴于各个人的舌头照片各异，颜色形状不尽相同，各种方法分割并不能适用于各种不同人的舌头。

4)暂时还没有用算法开发出应用程序直接用于临床使用。

发明内容

针对上述已知技术出现的缺陷和医生的需要，本发明的目的就是提供一种实用性较强的算法和基于该算法开发的程序，实现自动分割提取舌体的过程，而且使得到的结果精确，时间花费更少。同时算法的结果为用计算机进一步分析舌头的图像提供有效数据。

本发明的特征在于依次含有以下步骤：

步骤(1)，向计算机输入用数码相机采集好的舌头的彩色图片，其中，以舌体为主；

步骤(2)，在计算中设定以下两个调整参数：

分割速度，从慢到快分为9档，相应的取值范围为2到10的整数，代表降采样的倍数；

舌体轮廓曲线光滑度，从粗糙到光滑分为6档，相应的取值范围为1到5的整数，代表该曲线迭代的次数；

步骤(3)，利用分水岭变换方法按以下步骤对输入彩色图像的红色通道进行处理，得到粗略的舌体轮廓图：

步骤(3.1)，把所述图像中的每个像素点的红色分量取出，按照像素的排列方式组成一张灰度图表示的该图像的红色通道R(x，y)；

步骤(3.2)，按照设定的降采样倍数M对图像进行降采样，得到降采样后对应的表示在灰度图中的红色通道的值f(x，y)＝R(Mx,My)，其中，M为2-10之间的整数，(x，y)为降采样后的像素点的坐标，R(Mx,My)为降采样前图像在(Mx,My)点的红色分量；

步骤(3.3)，对所述灰度图f(x，y)，按下述利用标记控制的分水岭变换方法求出舌体的粗糙耸廓，建立舌体轮廓模型：

步骤(3.3.1)，计算灰度图f(x，y)的梯度图g(x，y)∶g(x，y)＝f(x，y)

步骤(3.3.2)，按以下步骤计算标记函数，用以标记灰度图f(x，y)中的目标和对应的背景：

步骤(3.3.2.1)，从图像RGB值计算饱和度St(x，y)：

St (x, y) = \frac{Max (R (x, y), G (x, y), B (x, y)) - Min (R (x, y), G (x, y), B (x, y))}{Max (R (x, y), G (x, y), B (x, y))}

其中，Max(R(x，y)，G(x，y)，B(x，y))≠0，St(x，y)为像素(x，y)点的饱和度，R(x，y)，G(x，y)，B(x，y)分别表示像素(x，y)点的红色、绿色、蓝色分量；

步骤(3.3.2.2)，根据设定的阈值St_th，对于每个像素点(x，y)的饱和度进行判别，得到标记图像f_m(x，y)：

饱和度阈值St_th＝0～1，取St_th＝0.4；

步骤(3.3.3)，对梯度图g(x，y)用标记图f_m(x，y)进行强制最小值，得到待分割图f_min(x，y)：

f_min(x，y)＝g(x，y)*f_m(x，y)，

f_m(x，y)为标记时，f_min(x，y)＝0，f_m(x，y)非标记时，f_min(x，y)＝g(x，y)；

步骤(3.3.4)，按以下步骤对灰度图f_min(x，y)用分水岭变换的方法从中提取出一片连通区域构成的舌体区域：

步骤(3.3.4.1)，在X-Y平面上，对灰度图f_min(x，y)沿着Y轴按给定的ΔY值作纵剖面图；

步骤(3.3.4.2)，在所属的纵剖面上，得到一条的采样点为横坐标，相应的灰度值作为纵坐标的曲线：

步骤(3.3.4.3)，在所属采样点一灰度图f_min(x，y)曲线上取出最小灰度值作为灰度阈值N，把f_min(x，y)分为大于灰度阈值和小雨灰度阈值上、下两个部分，构成的二值图像：

该二值图象T_N(x，y)的前一时刻的二值图像为

其中，Step为增加的步长，在1～5间取值；

步骤(3.3.4.4)，把步骤(3.3.4.3)所得两个相邻时刻的二值图象做异或操作，判断T_N-Step(x，y)≠T_N(x，y)时，S(x，y)是否等于1，T_N-Step(x，y)＝T_N(x，y)时，S(x，y)是否等于0，如果否定，则继续增加步长，再作判断，一直到满足为止，此时，S_j(x，y)是小分水岭，像素点(x，y)是分水岭的位置，像素点(x，y)对应的灰度值属于一片量筒的舌体区域的一部分；

步骤(3.3.4.5)，重复步骤(3.3.4.3)和步骤(3.3.4.4)，得到所有的分水岭，构成一个分水岭的集合S(x，y)＝S₁(x，y)·S₂(x，y)·...·S_K(x，y)；

步骤(3.3.4.6)，把位于图像中部的一片连通区域作为舌体区域，取出其边界V：

V＝{V₁，V₂…V_L}，V_i＝(x_i，y_i)，i＝{1，2，…L)，

其中，S(x_i，y_i)＝1，V上的点是分水岭变换计算出的封闭边界；

步骤(4)，把步骤(3)得到的边界V作为二维形变轮廓模型的初始曲线，按以下步骤计算出光滑的舌体轮廓曲线：

步骤(4.1)，对每个处于像素V_i周围设定领域中的像素V_i′，定义该像素点V_i′的能量项为：

E_i(V_i′)＝αE_int(V_i′)+βE_ext(V_i′)，

其中，E_int是依赖于轮廓形状的连续能量函数，使轮廓趋于光滑：

E_{int} ({V_{i}}^{'}) = \frac{1}{I (V)} {| | {V_{i}}^{'} - γ (V_{i - 1} + V_{i + 1}) | |}^{2},

I(V)是曲线V中各点的平均距离：

I (V) = \frac{1}{L} Σ_{i = 1}^{L} {| | V_{i + 1} - V_{i} | |}^{2},

对于封闭的曲线，

γ = \frac{1}{2 \cos (2 π / Len)},

Len表示该曲线的长度，

所属设定邻域，取值为5×5，

E_ext是图像灰度能量的函数，将变形轮廓向高灰度移动，等于向原图像的边缘移动，使E_ext(V_i′)＝-g(V_i′)，

α、β是加权常数，取2α＝β＝2，

若：E_i(V_i′)＜E_i(V_i)，则取V_i＝V_i′；

步骤(4.2)，根据设定的曲线光滑度参数，在设定的邻域内各个像素逐次迭代，使初始的舌体轮廓曲线收缩到使得能量E_i的位置，该收缩后最终得到的曲线内部的连通区域即舌体区域；

步骤(5)，把最终得到的曲线用所选降采样倍数扩充后，放入原彩色图像中，得到原彩色图像中的舌体轮廓曲线V＝{V₁，V₂…V_L}，V_i＝(M*x_i，M*y_i)，i＝{1，2，…L}，再把非舌体部分去处，填充白色，仅保存舌体部分，输出结果。

本发明采用了两步进行舌体的分割提取，第一步先得到舌体粗略的轮廓图，第二步再对轮廓线进行调整，得到精确的轮廓图。这种方法既保证了舌体区域正确性又保证了舌体轮廓曲线的精确性。该方法可以进行全自动的分割，不需要医生进行手动选择区域等复杂操作；另一方面，专门针对该算法开发了应用程序，应用程序可以和医生交互，由医生控制分割需要的时间和舌体轮廓曲线的光滑度。

经过对多幅医院采集的舌头图像的测试，该方法分割提取的图像结果非常准确。测试的图片的采样对象年龄层次广泛，同时也包括健康和疾病人群，充分说明该方法能够适应图像样本的多样性，没有受到个体差异的影响(图6，图7，图8分别为老人，小孩，青年照片和结果)。该算法和程序已经应用于临床，方法操作简单方便，结果稳定准确。

附图说明

图1：算法的流程图；

图2：标记控制的分割算法流程图；

图3：分水岭变换示意图；

图4：分水岭变换详细步骤示意图：

(A)阈值为N时的分割图像示意，

(B)阈值为N-Step时的分割图像示意，

(C)得到的小分水岭图像示意；

图5：二维形变轮廓模型算法示意图；

图6：该算法用于老人的舌头照片取得的结果；

图7：该算法用于小孩的舌头照片取得的结果；

图8：该算法用于青年的舌头照片取得的结果；

图9：取红色通道，降采样后的图像；

图10：梯度图像；

图11：饱和度图像；

图12：标记函数；

图13：待分割图像；

图14：分水岭位置；

图15：二维形变模型的初始轮廓；

图16：二维形变模型迭代50次后的效果图。

具体实施方式

该算法核心由Matlab实现。软件界面用Delphi编写。输入的舌头照片由中医大夫用数码相机拍摄，拍摄时病人平伸出舌头，方便医生将整个舌头摄入照片。

程序从数码相机中读入的彩色原始照片，大小为1024*960像素。得到红色通道图像后，取M＝5，即横、纵坐标都降采样5倍，得到如图9所示的图像。

求它的梯度图像，得到梯度图像如图10所示，舌头的边缘比起图9来被明显增强了。

求得其饱和度通道图像，如图11所示，舌头部分的颜色较暗，表示其饱和度较低。

将图11处理后，能够得到二值的标记函数，如图12所示，其中，白色部分为0，表示标记的部分，黑色部分为1，表示不是标记的部分。

用强制最小值将标记函数盖在梯度图上，得到待分割的图像，如图13所示。图13与图11相比，舌头部分的细节全部都被屏蔽掉，能够有效的避免分水岭变换的过分割问题。

图14为分水岭变换的结果，白色部分为0，黑色部分为1，黑色曲线即算法求出的分水岭的位置。

将图像中间的一块连通区域取出，按照顺时针得到区域的轮廓，即取得的二维轮廓形变的初始值。图15表示的是该轮廓叠加在图9上的效果图。

再对图15中的曲线用二维形变模型的方法进行迭代，迭代100次后，得到效果图如图16所示。与图15相比起来，曲线更加光滑，而且舌体的轮廓边缘更加准确。

最后将原图上的舌头分割出，其他部分用白色填充，方法完成。

本发明的具体实施步骤如下：

1)输入用数码相机采集好的舌头的彩色图片。

彩色照片需要以舌体为主，专用于舌诊的研究。

2)调整参数，选择分割需要的时间和舌体轮廓的光滑度。

可供调整的参数有两个：一个是分割速度，从慢到快分为9档；一个是曲线光滑度，从粗糙到光滑分为6档。第一个分割速度的参数用于第一步的分割(即步骤3)；第二个曲线光滑度的参数用于控制第二步的分割(即步骤4)。具体应用时也可以不进行该项操作，选择默认值即可。

第一个分割速度参数的取值范围是从2到10的整数，代表的意思是将图像降采样的倍数，将采样的倍数越大，处理的图像就会越小，从而处理的时间越快；第二个曲线光滑度的参数的取值范围是从1到5的整数，曲线迭代的次数。

3)利用分水岭变换对图像的红色通道进行处理，得到粗略的舌体轮廓图。

在计算机中，一幅彩色图像的每个像素点的颜色由三个字节来表示，三个字节分别代表该像素点的红色，绿色，蓝色分量。将每个像素点的红色分量取出，按照像素点的排列方式组成一张灰度图，即该图的红色通道R(x，y)。

按照步骤2中选择的降采样倍数对图像进行降采样，设确定的图像降采样倍数为M，按照步骤2的设定，M为2～10之间的整数。

f(x，y) ＝R(Mx,My)(x，y)为降采样后的图像点的坐标，R(Mx,My)为降采样前图像在(Mx,My)点的红色分量，f(x，y)为降采样后取出的红色通道的值。之后的处理都在该灰度图f(x，y)上进行。对该灰度图f(x，y)，利用标记控制的分水岭变换求出舌体的大致轮廓，建立舌体轮廓模型，流程图如图2所示，具体方法如下：

首先，取得f(x，y)的梯度图g(x，y)，用于增强图像的边界。

g(x，y)＝f(x，y)

然后，借助特征检测确定标记函数，标记函数的作用是标记出目标和对应的背景。根据拍摄图像中舌体部分的饱和度较低的特点，标记函数从图像的饱和度图中求出。从图像RGB值中得到饱和度的表达式为：

Max(R(x，y)，G(x，y)，B(x，y))≠0

St (x, y) = \frac{Max (R (x, y), G (x, y), B (x, y)) - Min (R (x, y), G (x, y), B (x, y))}{Max (R (x, y), G (x, y), B (x, y))}

(x，y)为点的坐标，St(x，y)为(x，y)点的饱和度，R(x，y)，G(x，y)，B(x，y)分别表示(x，y)点的红色、绿色、蓝色分量。

选取亮度阈值St_th，对于每个像素(x，y)，标记图像f_m(x，y)的定义为：

St_th的取值范围为0～1，根据经验，取St_th＝0.4。

对g(x，y)用图像f_m(x，y)进行强制最小值，得到待分割图像无f_min(x，y)。

f_min(x，y)＝g(x，y)*f_m(x，y)，即f_m(x，y)为标记时，f_min(x，y)为0， f_m(x，y)不是标记时，f_min(x，y)为本来梯度g(x，y)的值。

最后对f_min(x，y)采用分水岭变换。

借助涨水法来介绍分水岭技术。利用图3来进行讨论，为了简便，仅画出各个目标的一维剖面图。横坐标代表图像中的采样点，纵坐标代表该点的灰度值。假设有水从谷底涌出并且水位逐渐增高。如果从两个谷底涌出的水水位高过了其间的山峰，这些水就会汇合。如果要阻止水的汇合，就需要在山峰上修建水坝，且水坝的高度随着水位的上升而增高。这个过程随着全部山峰都被水淹没而结束。修建的所有水坝把土地分割成很多区域，这些水坝就构成了这片土地的分水岭。以图3而言，图中出现的三个波为全局，表示为从图像中任意取出的一个部分；每一个波为局部，在每个局部中包含一个波峰和一个波谷。

以上是分水岭的步骤，即以Step的步幅逐渐增加灰度阈值N，Step表示每次灰度增加的值，取值可选为1～5。每当N大于一个局部极大值，就把当时的二值图像(只区分陆地和水域，即大于灰度阈值和小于灰度阈值两部分)与前一时刻(即灰度阈值上一个值的时刻)的二值图像进行逻辑异或操作。

用M₁，M₂，…M_R表示f_min(x，y)中各个局部极小值，即谷底孔的像素位置，如图4(B)中所示，取灰度阈值N，得到二值图像

T_{N} (x, y) = \{\begin{matrix} 0 & f_{\min} (x, y) < N \\ 1 & f_{\min} (x, y) &GreaterEqual; N \end{matrix}

N-Step为前一时刻的灰度阈值，如图4(A)所示，此时二值图像为

T_{N - Step} (x, y) = \{\begin{matrix} 0 & f_{\min} (x, y) < N - Step \\ 1 & f_{\min} (x, y) &GreaterEqual; N - Step \end{matrix}

如果T_N-Step(x，y)内的两个连通区域，如图4(A)中所示的包含M_j和M_j+1的两个局部最小值的连通区域C_n(M_j)和C_n(M_j+1)在T_N(x，y)中合并成了一个连通区域，如图4(B)所示，将T_N-Step(x，y)和T_N(x，y)两幅图像进行异或操作，定义为：

S_j(x，y)＝T_N-Step(x，y)T_N(x，y)

如果T_N-Step(x，y)＝T_N(x，y)，则S_j(x，y)＝0；

如果T_N-Step(x，y)≠T_N(x，y)，则S_j(x，y)＝1。

(x，y)是点的坐标，S_j(x，y)标志的是分水岭，如果(x，y)是分水岭的位置，则S_j(x，y)＝1。如图4(C)中所示。

N的取值从波谷开始，到波峰结束，每次增加Step，就可以得到所有的分水岭为S₁(x，y)，S₂(x，y)…S_K(x，y)

分水岭变换的结果为所有分水岭的集合：

S(x，y)＝S₁(x，y)·S₂(x，y)·…·S_K(x，y)

得到的分水岭结果后，根据照片拍摄特点，一般舌头的位置都位于照片正中，取得位于图片中心的一片连通区域为需提取的舌体区域。取出这个区域的边界V。

V＝{V₁，V₂…V_L}，V_i＝(x_i，y_i)，i＝{1，2，…L}，其中S(x_i，y_i)＝1，即V上的点肯定是分水岭变换计算出的边界。

这个边界便是分水岭算法得到的舌体的边界曲线，也是二维形变轮廓模型算法的初始曲线。

4)用上述方法的轮廓作为初始用二维形变轮廓模型的迭代方法提取精确的舌体轮廓曲线。

分水岭变换得到的结果V＝{V₁，V₂…V_L}，V_i＝(x_i，y_i)，i＝{1，2，…L}为图像上一组排序的点的集合。处在轮廓上的点，通过解一个最小能量问题来迭代地逼近目标的边界。对每个处于V_i周围5×5的邻域中内的点V_i′，定义该能量项为：E_i(V_i′)＝αE_int(V_i′)+βE_ext(V_i′)

其中，E_int是依赖于轮廓形状的能量函数，选取其为连续能量，迫使封闭曲线变成圆环，使得轮廓更趋于光滑。

E_{int} ({V_{i}}^{'}) = \frac{1}{I (V)} {| | {V_{i}}^{'} - γ (V_{i - 1} + V_{i + 1}) | |}^{2}

在式中I(V)是归一化因子，是曲线V中各点的平均距离：

I (V) = \frac{1}{L} Σ_{i = 1}^{L} {| | V_{i + 1} - V_{i} | |}^{2};

对于封闭的曲线，取

γ = \frac{1}{2 \cos (2 π / Len)},

其中Len表示曲线的长度。

E_ext是依赖于图像性质的能量函数，选取其为图像灰度能量，将变形轮廓向高灰度移动。E_ext(V_i′)＝-g(V_i′)。实际上，由于我们选取的是梯度图像，变形轮廓应该是向原图像的边缘移动。

α、β是加权常数，在本算法中，考虑到舌头附近边缘比较明显，取2α＝β＝2，即外力作用为主导。

如果E_i(V_i′)＜E_i(V_i)，则V_i＝V_i′。

迭代的次数又步骤2中的第二个参数控制，参数选择为1～5，即迭代次数为10次到50次。迭代后曲线收缩到使得能量E_i最小的位置，即本算法的最后结果。该曲线内部的连通区域即我们需要提取的舌体区域。

5)将最终得到的曲线V＝{V₁，V₂…V_L}，V_i＝(x_i，y_i)，i＝{1，2，…L}扩充后放入到原本的彩色图像中，即V＝{V₁，V₂…V_L}，V_i＝(M*x_i，M*y_i)，i＝{1，2，…L}。将非舌体部分去除，填充白色，仅保存舌体部分，输出结果。

结果可供下一步的图像处理使用，也可供医生查看。

Claims

1.中医舌诊用的舌头彩色数码照片的舌体分割提取方法，其特征在于，依次含有以下步骤：

步骤(2)，在计算中设定以下两个调整参数：

步骤(3.2)，按照设定的降采样倍数M对图像进行降采样，得到降采样后对应的表示在灰度图中的红色通道的值f(x，y)＝R(Mx，My)，其中，M为2～10之间的整数，(x，y)为降采样后的像素点的坐标，R(Mx，My)为降采样前图像在(Mx，My)点的红色分量；

步骤(3.3)，对所述灰度图f(x，y)，按下述利用标记控制的分水岭变换方法求出舌体的粗糙轮廓，建立舌体轮廓模型：

步骤(3.3.1)，计算灰度图f(x，y)的梯度图g(x，y)：g(x，y)＝f(x，y)

步骤(3.3.2.1)，从图像RGB值计算饱和度St(x，y)：

St (x, y) = \frac{Max (R (x, y), G (x, y), B (x, y)) - Min (R (x, y), G (x, y), B (x, y))}{Max (R (x, y), G (x, y), B (x, y))}

饱和度阈值St_th＝0～1，取ST_th＝0.4；