阈值确定方法

一、问题重述

图形（或图像）在计算机里主要有两种存储和表示方法。矢量图是使用点、直线或多边形等基于数学方程的几何对象来描述图形，位图则使用像素来描述图像。一般来说，照片等相对杂乱的图像使用位图格式较为合适，矢量图则多用于工程制图、标志、字体等场合。矢量图可以任意放缩，图形不会有任何改变。而位图一旦放大后会产生较为明显的模糊，线条也会出现锯齿边缘等现象。

矢量图从本质上只是使用曲线方程对图形进行的精确描述，在以像素为基本显示单元的显示器或打印机上是无法直接表现的。将矢量图转换成以像素点阵来表示的信息，再加以显示或打印，这个过程称之为栅格化（Rasterization），见图1。

栅格化的逆过程相对比较困难。假设有一个形状较为简单的图标，保存成一定分辨率的位图文件。我们希望将其矢量化，请你建立合理的数学模型，尽量准确地提取出图案的边界线条，并将其用方程表示出来。

二、问题分析

本题的要求是完成位图的矢量化，通过建立合理的数学模型，将一个有一定分辨率的位图文件尽量准确地提取出图案的边界线条，最终将位图用方程的形式表示出来。解决本问题的流程图见下图。首先，通过MATLAB读取位图的各个像素的像素值（0-1），得到位图各个点的灰度值，通过最大类间方差法和最大熵法确定阈值，完成灰度的二值化，使各个像素点的灰度值全部由0或1表示。其次，将位图的轮廓通过合适的算法提取出来，根据特征值对轮廓进行拟合。最后，根据拟合的函数完成位图的矢量图，完成其矢量化过程，并通过对比矢量图和原始位图对应的。

三、问题假设及符号说明

3.1问题假设

3.2符号说明

四、模型建立

4.1模型准备

本题要求将一个形状较为简单的图标，保存成一定分辨率的位图文件，即将位图矢量化。

阈值：指释放一个行为反应所需要的最小刺激强度，本文指像素点灰度值二值化的临界值。 4.2阈值的确定方法 4.2.1最大类间方差法

最大类间方差法的基本思想是将待分割图像看作是由两类组成的整体，一类是背景，一类是目标[6]。用各个像素点像素值的方差来衡量目标和背景之间额差别。方差越大，即两像素值的区别越明显，不能分在同一类别中。本问中，为使得目标和背景两类的类间方差最大的灰度级别即认为是最佳阈值。一副大小为

N M ?个像素、灰度级别为L 的图像，若设图像中灰度级为i 的像素个数为i N ，则灰度级i 的概率为

M

N N P i

i ?=

根据最大类间方差法的分割思想，最佳阈值公式为

()()[]

2

02

01

0*max arg w w P w w P T b b a a L t -+-=-≤≤

其中，

∑==T

i i a P P 1

∑+==

L

T i i

b P P 1

∑==T

i a i

a P P i

w 1 ∑+==

L

T i b

i

b P P i

w 1

∑==L

i i iP w 1

一维最大类间方差法是建立在一维灰度值分布图基础上的确定阈值的方法，而一维灰度值分布图没有考虑图像的空间信息，仅仅表示出了各个像素点独立的灰度分布，但是图像中的每个像素与其相邻的像素是有相关性的，所以在以为灰度值分布图基础之上的阈值分割方法所建立的分割模型考虑的因素不全面。因此，在对图像进行分割时，容易受到噪声的干扰，得到的分割结果和阈值结果不精确，将影响后续位图矢量化的准确度。

为了增强算法的抗噪性能、提高其确定阈值结果的精确度，本文将最大类间方差法从一维算法推广至二维算法。二维最大类间方差法是建立在二维灰度值分布图的基础上计算阈值的算法，下面对其进行介绍。

设一幅大小为N M ?个像素的图像灰度级为L ，()y x f ,为图像在()y x ,的灰度值，()y x g ,是以()y x ,为中心，k k ?邻域内的平均灰度值，其表达式为

()()∑∑=-==-

=++=

222

2

2

,1,k m k m k n k n n y m x f k y x g

式中，k 为邻域的大小；

()y x g ,为邻域内的平均灰度值；

()n y m x f ++,为图像在点()n y m x ++,的灰度值； M 为图像的宽度； N 为图像的高度。

()y x f ,和()y x g ,即形成了一个二元组()()()y x g y x f ,,,，令ij N 为图像中

()i y x f =,且()j y x g =,的像素的个数()1,0-≤≤L j i ，则其联合概率密度为

N

M N P ij ij ?=

式中，ij P 为图像中()i y x f =,且()j y x g =,的像素联合概率密度； ij N 为图像中()i y x f =,且()j y x g =,的像素的个数。

由最大类间方差法求出的阈值向量()t s ,会将直方图分为四个部分，分别为目标部分、背景部分、边界部分、噪声部分。在大多数情况下，近似的认为边界部

分和噪声部分的概率为0，因此这种近似在多数情况下是合理的。

图像目标出现的概率为

∑∑===s i t

j ij P P 11

图像背景出现的概率为

∑∑+=+==

L s i L

t j ij

P

P 11

1

式中，ij P 为图像中()i y x f =,且()j y x g =,的像素联合概率密度； L 为图像灰度级。

图像总均值矢量为

()'????

??='=∑∑∑∑====L i L j L i L j ij ij T T T jP iP 111110,,μμμ

图像目标和背景额均值矢量为

()'???

?

??='=∑∑∑∑====s i t j s i t j ij ij jP iP 111101000,,μμμ

()'???

?

??='=∑∑∑∑+=+=+=+=L s i L t j L s i L t j ij ij jP iP 111111101,,μμμ

因为假设边界区域和噪声区域的概率为0，因此有

110≈+p p

1100μμμp p T +=

式中， 0μ和1μ分别为目标区域和背景区域的均值矢量。 类间方差定义为

()()()()??????'

--+??????'--=T T r r B p p S μμμμμμμμ111000

取类间方差的迹作为测度，有

()()()

002

112001p p p p trS T j T i B --+-=μμμμ

式中，

∑∑===s i t

j ij i iP 11

μ

∑∑+=+==

L

s i L

t j ij

j jP 11

μ

式中，ij P 为图像中()i y x f =,且()j y x g =,的像素联合概率密度；

L 为图像灰度级。 最佳阈值公式为

()()t s trS t s B

L

t L s ,max arg ,11*

*

≤≤≤≤=

式中，B trS 为测度。

遍历图像的灰度级，求使目标和背景的类间方差的迹最大的灰度级，即得最佳分割阈值()**,t s 。 4.2.2最大熵法

由熵的定义可知，当一个信息源中所有的事件以相同的概率出现时，这时候信息熵大，因为信息源的不确定性大，所以信息熵大。用最大熵法求最佳阈值，就是求一个分割阈值使得目标和背景两类的信息熵之和最大，一幅大小为

N M ?个像素、灰度级为L 的图像，设图像中灰度级为i 的像素个数为i N ，则灰度级i 的概率为

M

N N P i

i ?=

图像目标区域的熵和图像背景区域的熵可以表示为

∑

=-=T

i t

i t i f P P P P H 1log ∑

+=---=L

T i t i t

i b P P

P P H 11log 1 式中，f H 为图像目标区域的熵； b H 为图像北京区域的熵； ∑==T

i i t P P 1

最大熵法的最佳阈值公式为

()()[]

T H T H T b f +=max arg *

遍历图像的灰度级，求目标区域的熵和背景区域的熵的和最大的灰度级，即得最佳分割阈值*T 。

4.2.4 确定阈值方法总结

将以上提出的累积剩余熵阈值分割算法用于实际图像分割，将分割结果与经典的最大类间方差法和最大熵法进行了比较。实验中用了人、直升机、车、帆船、船和发电站图。表1是对分割结果的主观评价，表2是三种阈值分割方法的最佳分割阈值的比较。

表1 分割结果主观评价表

图像 最大类间方差法

最大熵法 最大累积剩余熵法

人 差 好 好 直升机 差 好 好 车 差 好 好 帆船 一般 一般 一般 船 差 一般 好 发电站

差

一般

好

表2 最佳分割阈值比较表

图像 最大类间方差法

最大熵法 最大累积剩余熵法

人 129 165 153 直升机 82 94 94 车 65 85 77 帆船 134 156 137 船

99

126

144

发电站 100 177 165

由上表分析可知，图像的不同对应的最佳分割阈值的确定方法不同，如：当图像为人、直升机或者汽车时，最大熵和最大累积剩余熵确定阈值的方法比最大类间方差法结果更精确。 4.3二值化

根据像素点的灰度值()y x f ,和阈值T ，可将图像二值化。

()()()??

?<≥=T

y x f T y x f y x f ,0

,1

,

式中，()y x f ,为像素的灰度值，T 为最佳分割阈值。 4.4基于细化算法的轮廓提取 4.4.1基础术语

为简便的描述图像中像素之间的相互关系，本文引用了邻域和连通的概念[8]。

●4邻域和4相邻

对于图像中的某个像素点()n m p ,，则与之在水平方向（左和右）和垂直方

向（上和下）相邻的4个像素点坐标分别为()n m ,1-，

()n m ,1+，()1,-n m ，()1,+n m ，则这4个像素点组成了像素点()n m p ,的4邻域，表示为()p N 4。而这4个像素点在位置上就与像素点()n m p ,相邻。

● ●

○

● ●

图 4邻示意图 图 4邻坐标关系图

●8邻域和8相邻

若取像素()n m p ,四周的8个像素点作为相邻点，则像素点()n m p ,的这8个

相邻点就构成了8邻域，用()p N 8表示。 ● ● ● ● ○ ● ●

●

●

图 8邻示意图 图 8邻坐标关系图

●像素间的邻接

像素的相邻仅说明了两个像素在位置上的关系，若再加上取值相同或相近，则称两个像素邻接。

①两个像素点p 和q 邻接的条件

像素点()n m p ,和()t s q ,位置上满足相邻，即 4相邻：()()q N n m 4,∈或者()()p N t s 4,∈ 8相邻：()()q N n m 8,∈或者()()p N t s 8,∈

且像素点()n m p ,和()t s q ,灰度值相近，即称为灰度值相近（似）准则，即

V p ∈ 和V q ∈，其中{},...,21v v V =

②邻接的类别

4邻接： V q p ∈,且()()p N t s 4,∈，则称像素点()n m p ,和()t s q ,4邻接，记为

q p ?→?4；

8邻接： V q p ∈,且()()p N t s 8,∈，则称像素点()n m p ,和()t s q ,8邻接，记为

q p ?→?8。

两种邻接的关系：4邻接必8邻接，反之不一定成立。

两种邻接及其关系见下图，相似性准则为{

}1=V ，其中像素点p 和q 4邻接，也

8邻接，点q 和r 8邻接但非4邻接。

图 邻接关系图

●像素间的连通 ①通路

设()n m p ,与()t s q ,之间的各像素点形成的连线L 为

()()()()()()()(){}t s n m n m n m n m n m n m q p L n n i i i i ,,,...,,,,,...,,,,,.111100===--

若()11,--i i n m 与()i i n m ,邻接()n i ≤≤1，则()q p L ,称为()n m p ,与()t s q ,之间的一条通路，N 为通路长度。与链接一样，通路也分为4通路和8通路。 ②连通性

若S 是图像中的一个子集，S q p ∈,，且存在一条由S 中像素组成的从p 到q 的通路，则称p 在图像集S 中与q 连通，连通也分为4连通和8连通。

图连通示意图

4.4.2 二值图像的细化

所谓细化，就是在保持原图像拓扑结构的情况下，将线性图转成单像素宽的骨架。细化算法是一种特殊的多次收缩的过程，由于它在收缩过程中不会消去只有一个临接点的边界点，因此可以保证收缩过程中的连通性，细化后的曲线位于区域中央，一般采取由左至右、由上至下的细化过程，反复多次地消去具有两个

以上临接点的边界点。其基本思想是对某一像素p，观察其邻点像素的分布情况，在不破坏图形连通性的情况下，将p点的值由1变为0。下图为细化过程实例。

(a) 原始图像(b) 消去上边界后的区域图像

(c)消去下边界后的区域图像(d)消去左边界后的区域图像(e)消去右边界后的区域图像

图细化过程示例图

4.5拟合曲线

通过对图像的细化得到一个单像素宽的骨架轮廓，然而，若将组成骨架的一系列点全部拟合到一个函数中，拟合结果的精确次数将会过高，不便于计算机的运算。因此，为简化运算、降低函数的次幂，我们对构成单像素宽骨架的所有点进行了基于特征的分类和提取，为保证结果的简洁性，根据不同类型点的特征对

骨架进行分割，分别拟合各段轮廓，得到幂级数较低的拟合函数曲线。 4.5.1提取分类点

根据构成单像素宽骨架的点的特征，将点分为三类：间断点、分割点、特征点。 ●间断点

间断点为轮廓中的突然出现的尖点，类似于数学概念中的不可导点，由于图像中的像素点都是离散的，因此对此类点也进行离散化处理。本文引用差分的方法，用一阶差商来近似一阶偏导数，一阶差商有三种形式：一阶前差商、一阶后差、一阶中心差商，分别记为

x p p x p i

i ?-=

??+1 x p p x p i i ?-=

??-1

●分割点

分割点为轮廓中曲率突变的点，类似于数学概念中刻画离散点的曲率，可用于将骨架轮廓分割成段。曲线中变化剧烈的点的曲率较大，在保证精度的前提下，变化曲率较大的地方要特殊保留，可以通过估算曲线上各点的余弦函数值来间接地估算该点的曲率，以此得到曲线上变化剧烈的点作为明显分割点进行先行保留。

设i p 点是轮廓曲线连通段上点序中第i 个离散点，k i p -和k i p +分别是距离i

p 点步长为k 的前后两个点，则i p 点的曲率值可用下列公式估算。

()()()()()()()()2

222cos i k i k i i i k i k i i i k i k i i i k i k i i i y y y y x x x x y y y y x x x x -+--+---+--=

+-+-+-+-θ

根据阈值大小，将曲率值大于阈值的点作为分割点。 ●特征点

特征点为轮廓中能最大限度代表整个轮廓形状的点，类似于近似轮廓组成点。为避免拟合函数中出现有高度波动性的高次多项式，特征点的选取可简化轮廓，最终降低拟合函数的次数。

以图为例进行拟合，特征点的选取要经历如下步骤。

确定一个拟合误差ε，连接始点A 和终点H ，找出距离AH 直线的最远点C ，若C 到AH 的距离小于ε，则AH 即为满足要求的拟合矢量，若C 到AH 的距离大于ε，则C 到AH 之间的一个锚点，可用AC 和CH 两个矢量代替AH 矢量。用同样的方法依次分别对AC 和CH 处理，直到所有离散点与其拟合矢量的距离均小于预定的误差ε为止，最终以上满足预定误差ε条件的离散点即为单像素宽

轮廓的特征点。

图基于距离判断的特征点提取示例

用有限次端点迭代的方法，找出单像素宽轮廓图中满足一定精度的拟合矢量。以图为例，具体的取特征点并矢量化的过程如下[7]：（16

≤ε）

1≤

①以AH为拟合矢量，找出距离AH的最远点C，C距AH>ε，则C为锚点；

②以AC、CH代替AH，先找出AC区段内所有点位距AC的最远点位B，B距AC>ε，则B仍为锚点；

③以AB、BC代替AC，先找出AB区段内距AB最远点位并发现该点位距AH>ε。则B点为一个顶点。再找BC区段内距BC的最远点位并发现该点距BC<ε，则C点为顶点。至此AC区段已处理完毕；

④找出CH区段内距CH的最远点位D，D距CH>ε，则D为锚点；

⑤以CD、DH歹徒CH，先找CD区段内距CD的最远点位并发现该点距CD<ε，则D为顶点。再找DH区段内距DH最远点位E，E距DH>ε，则E为锚点；⑥以DE、EH代替DH，先找DE区段内距DE的最远距离点位并发现该点距DE<ε，则E为顶点。再找EH区段内距EH的最远距离点F，F距EH>ε，F为锚点；

⑦以EF、FH代替EH，先找EF区段内距EF的最远点位并发现该点距EF<ε，则F为顶点。再找FH区段内距FH最远点位G，G距FH>ε，则G为锚点；

⑧以FG、GH代替FH，先找GH区段内距FG最远点位，发现该点距FH<ε，则G为顶点，再找GH区段内距GH最远点位，发现该点距GH<ε，则H为最后一个顶点。

可见，算法每次迭代所傲的工作都是相同的：找出拟合矢量区段内的最远点位，求出该点到拟合矢量的距离，将这个距离与ε比较,然后分两种情况，这个距离若小于ε，则拟合矢量的末点为顶点，若大于ε，则为锚点，因此易于编程。在寻找拟合矢量区段内的最远点位时，仅对该区间内的所有点位与拟合矢量之间进行水平（或者垂直）方向的距离比较，不必求出每点到拟合矢量的真正距离。

4.5.2拟合函数

在得到特征点的基础上，我们建立了基于最小二乘法、Bezier 曲线和B 样条的曲线拟合模型。由于最小二乘法的原理较为通俗易懂，本文就贝奇尔曲线和B 样条的拟合方法作重点介绍。 ●Bezier 曲线拟合法

一条n 次Bezier 曲线[9]由n+1个控制点决定，依次连接这n+1个控制点即得到特征多边形，如下图所示。

由Bezier 曲线定义可得

()()∑==n i i i n V u J n r 0

,v v

式中，n 为Bezier 曲线次数；

i 为特征多边形顶点的序号，n i ≤≤0； u 为参数，10≤≤u ；

i V v

为特征多边形顶点的位置矢量；

i n J ,为伯恩斯坦基函数。

()()i

n i

i n i n u u C u J --=1,

其中，i

n C 为组合数，

()!

!!

i n i n C i n -=

为表述简便，可将以上公式表示成矩阵的形式，即

()[]

T n n n V V V Mn t t t n r ??

????=-v v v v ,...,,1,,...,,101 其中，=Mn 是n 阶方阵，元素ij b 为

()?

?

?+≤+-=----else

n j i C C b j i

n i

n n j i n ij 01

1

已知控制点的位置矢量是使用Bezier 曲线的前提，由于已知N+1个型值点分别为()N j P j ,...,2,1,0=，可用逆算法求解出控制点，即

(),0

n

n i

j

i

j

i J u V P ==∑

写成矩阵的形式即为

()()

()()()()()()??????

???????????????????

??

?------=--------n n n n n

n n n n n n n

n

n

n n n n n n n n n n V V V C C C

C

C C C C C C

C C t t t u r v v v v ...0...

1............0...111...111,,..,,10000011

12011111

1

该方程有n+1个未知数，N+1个方程，若n=N ，则方程有唯一解；若n>N ，则方程有无穷多解，即控制点可随意选取；若n

根据上述分析可知，特征点的个数一般小于4个，而且求得的Bezier 曲线要与两侧的曲线光滑拼接，故令n=N+2，将多余的两个方程用于边界条件的求解，从而使得方程有唯一解。

边界条件：设与()u r v

相邻的两条曲线如下，且[]1,0∈u

()()∑==m i i i m P u J u s 0

,v v

()()∑==l i k i l Q u J u w 0

,v v

①位置连续条件：

()()?????==0

1

10Q r P r v v v v ②斜率连续条件

()()()()

?????'=''='0101w r r s v v v v βα 其中，βα, 为参数。当1≠≠βα时，即切矢共线，但模长不同，可满足数学上可到的条件；当1==βα时，拼接时切矢量相等，即曲率相等，此时拟合效果最好。 ●B 样条拟合法

与Bezier 曲线拟合法类似，B 样条曲线方程为

()()∑==n i i k i V u N u r 0

,v v

式中，k i N ,为B 样条的基函数；

i 为基函数序号； k 为基函数次数。

()()()()

???

???

??--+--=????≤≤=?-++++++-++t N t t t t t N t t t t t N else t t t t N k i i k i k i k i i k i i k i i i i 1,11111,,1

,01 规定

00

=。其逆运算与Bezier 曲线类似， ()j i n

i i k i P V u N v

v =∑=0

,

同理，也多选两个控制点来补充端点条件。

五、模型求解及结果分析

5.1模型求解

●Step 1：将示例位图导入MATLAB ，并读取位图的像素灰度值；

●Step 2：在得到灰度值的基础上，选用合适的阈值确定方法，进行图像二值化；

●Step 3：将二值化后的图像进行细化，得到单像素宽的图像轮廓骨架；

●Step 4：将细化像素提取特征点及分割点；

●Step 5：用两个分割点的特征值进行Bezire曲线拟合，得到有两个分割点的函数方程；

●Step 6：将所得到的分割点函数方程进行拼接，完成矢量化。

5.2模型结果及分析

实例位图如下图所示。

图示例原始位图

本例在求解过程中，使用了最大类间方差法和最大熵法两种方法来确定最佳分割阈值，从而进行图像的二值化处理，两种方法得到的二值化图如下。

图最大类间方差法二值化图图最大熵法二值化图

两种方法细化后的图像如下图。

图最大类间方差法细化图图最大熵法细化图

为标记提取的间断点和分割点，我们将细化后的图进行颜色的黑白反转，使背景为白色，轮廓为黑色，以便更加明显地标记拟合需要的各类点。颜色反转后的图像如下。

图最大类间方差法细化颜色反转图图最大熵法细化颜色反转图

由上图对比可知，两种方法得到的二值化图和细化图与原始图像的差距都不大，且都较好的反映了原始图像的轮廓，但是，最大类间方差法的二值化图像和细化图像中的噪点比最大熵法的二值化图像和细化图像中的噪点多，因此，为更准确的进行特征点的提取和曲线拟合，我们采用最大熵法的图像进行后期处理。提取间断点和分割点后的图像如下。

图间断点分割点分布图

六、模型评价

6.1模型优点

1.

6.2模型缺点

1.

七、模型优化与推广

八、参考文献

[1] Pal S K ,Pal N R. Entropic thresholding [J].Signal Processing,1989,16(2):97-108.

[2] 潘喆,吴一全.二维指数熵图像阈值选取方法及其快速算法[J].计算机应用, 2007,27(4)：982-985.

[3] Sahoo P K, Wilkins C, Yeage J. Threshold selection using Renyi's

entropy[J].Pattern Recognition,1997, 30 ( 1) : 71-84.

[4] Abutaleb A S. Automatic Thresholding of Gray-Level Picture Using

Two-Dimensional Entropies [J].Computer Vision ,Graphics ,and Image Processing ,1989,47(1):22-32.

[5] Rao M ,Chen Y ,Vemuri B C,etal. Cumulative residual entropy:A new measure of information [J].IEEE https://www.360docs.net/doc/7b14745508.html,rm.Theory,2004,50(6):1220-1228.

[6] 郭臻, 陈远知. 图像阈值分割算法研究[J]. 中国传媒大学学报: 自然科学版, 2008, 15(2): 77-82.

[7] 宁飞. 计算机图象处理中点阵图矢量化的研究及其在图形识别中的应用[J]. 微型机与应用, 1988, 5: 004.

[8] 王晓丹, 吴崇明. 基于MATLAB 的系统分析与设计: 图像处理[M]. 西安电子科技大学出版社, 2000.

[9]计算机辅助几何造型技术[M]. 科学出版社, 2009.