根据子结点像素计算父结点像素时常用的是灰度平均法。我的应用中有所不同：如果四个像素中白色占了3块以上则父结点像素赋值为白色；如果黑白2：2对角分布也赋值为白色；如果黑色占了三块以上则赋值为黑色；如果黑白2：2各占据一侧则赋值为黑色。这样的好处是可以消除一部分不紧密的连接（2：2对角分布），而保留所有可能的紧密的连接（2：2同侧分布）。

分合算法的思路是^（5）：

1．任取一层数据结构，检验其均匀性H，若对于该结构中某一区域R有H(R)=false，则分裂这个区域为4个子区；若4个彼此相邻属于同一个父结点的区域R_k1，…，R_k4满足H(R_k1∪R_k2∪R_k3∪R_k4)=true，则合并它们成为单一区域。当不再有区域可分裂或合并时，停止分割。

2．如果二任意相邻区域R_i和R_j（可以不属于同一个父结点；可以有不同的大小）使得H(R_i∪R_j)=true，则合并它们。

 

2.2.1 分割算法

分割算法的描述如下：

1.把金字塔数据结构中起始某个中间层的数据块编码全部压入工作堆栈RgStack。

2.从RgStack弹出一个数据块编码Code。如果Code代表的数据块不均匀，就把Code分裂为四个数据块，再把颜色为黑色的那些块的编码压入堆栈RgStack。如果均匀，或者Code代表的已经是不可分的像素，那么把Code压入中间结果堆栈RgCode。

3.反复执行2，直至堆栈RgStack空。

有两点值得解释一下。首先，我在程序中采用整形数来保存数据块编码，这样较字符串而言更加方便了计算（子结点或父结点编码的计算，数据块坐标的计算等等）。其次，我使用了递归方法产生中间层所有数据块的编码,相对于多重循环而言，不仅使程序书写简洁，而且起始的中间层号可以方便地修改，而多重循环结构要修改起始层则必须加或减循环的层数。

 

2.2.2 合并算法

合并算法的描述如下：

1.把中间层里所有均匀的数据块的编码存入工作队列RgA。

2.从RgA出队一个编码Code。在RgA中寻找和Code属于同一个父结点的其他数据块。如果找不齐所有四个，那么把Code和找到的都从RgA删除，而压入中间接果数组RgCode。如果找齐了，那么把这四个编码都从RgA删除，而把父结点的数据块编码入队到RgA中。

3.反复执行2，直至队RgA空。

需要说明的一点：步骤1虽然是属于合并这部分的，但是在执行分裂时就可以方便地做掉，节省代码和时间。

 

2.2.3 相邻归并算法

相邻归并算法的描述如下：

1.从RgCode中弹出一个编码，压入工作堆栈RgStack。

2.从RgStack中弹出一个编码Code，将它入队到RgA。扫描RgCode，把所有与Code相邻的数据块的编码都压入堆栈RgStack。

3.反复执行2，直至栈RgStack空。此时RgA里面存有的一系列数据块编码就是属于同一个特征块的，可以用来计算这个特征块的形状属性并保存起来。具体计算方法见2.4“形状属性的计算”。

4.反复执行1-2-3，直至栈RgCode空。

这里有两点值得说明。首先，相邻归并算法是三个算法中最耗时的，占去了分合算法中绝大部分的运行时间（对一幅182×248单人脸灰度图像做试验，分裂，合并，相邻归并共用3224毫秒，其中相邻归并占了3074毫秒）。而且其耗时随问题规模是非线性增长的。对于小而简单的图像，分合算法可以很快（3秒钟），可是对于较大或较复杂的图像，分合算法往往耗时很长（20～40秒）。这么大的差别就是来自于非线性增长的计算量。其次，判断两个编码所代表的数据块是否相邻的函数在相邻归并算法中被频频调用，它的效率会大大影响整个算法执行的时间。我曾经采用计算两个数据块中心的x向和y向距离，然后与半边长之和比较的方法来判断数据块是否相邻，方法直观，代码较少，但是由于涉及了浮点运算，并且对每一个编码都要做完整的一套计算，导致运行速度很慢。现在采用的是分析两个编码在各位上的组合方式，从大尺度开始一位一位进行检查，只要发现这两个数据块在某一尺度上不可能相邻，就马上做出不相邻的判断，除非它们能够一直保留到被证实符合一种相邻的编码组合。具体方法描述较繁，这里省略。采用这种方法之后，虽然代码长度增加了很多，但是执行效率却大大提高了。

 

2.3 成组算法

引入成组算法出于如下事实：一些本属于同一个脸部器官的特征块在分合之后仍然是分开的。分开的特征块的中心、大小、取向等等参数相对于一个代表完整的器官的特征块而言都是错误的，导致后面的几何匹配得不到正确结果。成组算法就是用来把挨得足够近的可能本属于同一个器官的特征块合并起来（尽管它们并不是连通的）。算法实现参考了Shi-Hong Jeng, Hong Yuan Mark Liao等^（1）的工作。算法有两个特点：1. 合并半径取决于待合并的特征块的大小（一个重要参数－最大合并半径N_max事先指定），特征块越大，则合并半径越小；2. 如果待合并的两个数据块中任一个大于N_max­，那么它们必须位于同一主轴上（主轴的定义在形状属性的计算中介绍）才能合并。下面的伪代码描述了成组算法，其中N(i)代表第i个块B_i包含的像素个数，N_b代表图像中特征块的总数。

 

Procedure Grouping

for N := 1 to N_max do

  R := 8 – 6 * (N-1) / (N_max-1)

  for i :=1 to N_b do

  if N(i) = N then

    “找到所有位于B_i的r邻域内的特征块{B_j}”

  if {B_j} ≠φ then

    “从{B_j}中找出离B_i最近的块B_k”

  if  (N(k) < N_max) or (B_k与B_i位于同一根主轴上) then

    “把B_i和B_k合并”

  end for

end for

end Grouping

 

合并后的特征块的各个形状属性都要重新计算，具体方法见2.4“形状属性的计算”。

 

2.4 形状属性的计算

二值化后的原始图像直接用黑和白的像素表示。在用分合算法分割之后，每一个连通的黑色特征块就由一组数据块编码来表示，编码所代表的数据块（可能是来自金字塔数据结构中不同的层次）在一起组成了该特征块。但是这样仍然不能方便地进行几何匹配。为了进行几何匹配，需要知道每一个特征块的位置，大小，倾角及大概的形状。这些参数统称为形状属性。从组成一个特征块的数据块编码来计算这个特征块的形状属性是继分合算法之后的重要步骤。由于它与分合算法结合更紧密（因为不同的分割算法导致不同的计算形状参数的方法），所以把它放在分合算法这部分来介绍。

图像用f(x,y)来表示，x，y分别是横、纵坐标（已经离散化，所以下面都使用求和而不是积分），f(x,y)是该坐标处的灰度值。用字母B代表特征块。设特征块由一组正方形数据块{B_i}组成，B_i包含的像素个数为N_i，边长为D_i（N_i = D_i²），其中心（也就是重心）坐标为。N_i和_<![endif]>可以从数据块编码通过计算得到，方法较简单，这里就不赘述了。形状属性的计算就是要把特征块的形状属性都用N_i（或D_i）和表达出来。

1．特征块的位置：

特征块的位置用它的重心来表示。重心坐标计算公式如下:

在我的应用中，像素都是黑白的，即f(x,y)只有两种取值：0（黑）或者255（白），而特征块所包含的像素一定是黑色的（这是特征块的定义决定的）。所以在公式中完全可以把f(x,y)省略（为了方便，下面的公式中都直接略去了f(x,y)），直接写作：

经过简单推导得到用N_i和<表达的重心坐标：

，

 

其中是特征块包含的像素个数，即特征块的大小。

2．特征块的取向：

特征块的取向用特征块的主轴的倾角（规定取值范围）来表示。特征块的主轴定义为使转动惯量I(θ)取值最小的转轴⁽⁶⁾。可以证明（证明略），主轴一定穿过重心，而且使转动惯量取极值的倾角θ满足条件⁽⁴⁾：

_{，其中μ­11，μ20，μ02是三个二阶中心矩}

由此可以得到：

_<![endif]>，k为整数

由于θ取值范围的规定，k的所有可能取值只有-1，0，1。

算法中对k的所有三个可能取值分别计算θ，舍弃不在_<![endif]>内的值。然后分别检查转动惯量二阶导I^’’(θ)的符号，使得I^’’(θ)>0的才是能够令I(θ)取极小值的主轴倾角θ。检查I^’’(θ) 的符号所依据的公式是：

计算主轴倾角θ要用到中心矩。中心矩的一般定义和详细讨论可参见Robert M. Haralick 和 Linda G. Shapiro 的 Computer and Robot Vision⁽⁶⁾。

特征块B的（j+k）阶中心矩μ_jk为（注意这里同样略去了f(x,y)）：

_<![endif]>

经简单计算可以得到：

因此问题转化为求_<![endif]>。

这一部分的推导比较繁琐，结果如下：

有了这些公式，特征块的取向就可以由特征块所包含的数据块的编码计算出来了。

3．特征块的外接矩形

特征块的外接矩形定义为长边平行于主轴且与特征块外接的矩形。

为了求得外接矩形，对特征块做坐标变换：把坐标原点挪到重心处，再逆时针旋转坐标系θ，得到的新坐标系记为po^’q。在新坐标系下找出所有数据块顶点（每个数据块有四个顶点）坐标的最值p_min，p_max，q_min，q_max。易知在新坐标系下外接矩形的四个顶点坐标为（p_min，q_min），（p_min，q_max），（p_max，q_min），（p_max，q_max），特征块重心到矩形四条边的距离分别是L1 = - p_min，L2 =  p_max，W1 = - q_min，W2 =  q_max。将坐标系反变换回去（先顺时针旋转θ再平移），就得到了原坐标系下特征块的外接矩形的四个顶点坐标。参见图4。

最后讨论一下两个特征块（“物体”）合并时形状属性的计算。设B₁和B₂这两个特征块合并成了B。

显然有N = N₁ + N₂；_<![endif]><<![endif]>。

由刚才得到的公式：

_<<![endif]>

刚才i是对数据块计数的，现在推广到对特征块计数。可见如果把每个特征块的∑xy，∑x²，∑y²这三个数据保存在“物体”属性里，那么合并以后的特征块的∑xy，∑x²，∑y²很容易就能计算出来，其主轴倾角也就容易计算了。

合并后“物体”的外接矩形理论上应该由所有组成这个“物体”的数据块来计算。由于在“物体”属性里不能再保留所有成员数据块的编号（否则存贮量太大），所以合并后“物体”的外接矩形无法精确求出。可以把各个子“物体”的外接矩形的外接矩形近似当作合并后“物体”的外接矩形，如图5。计算的方法与单个特征块外接矩形大同小异，这里就不重复了。

TIFY: inter-ideograph; TEXT-INDENT: 25pt; TEXT-ALIGN: justify">分合完成，物体属性都计算完毕之后，就可以在图像中标出所有的特征块。在我的程序中，每个特征块用一个外接矩形和一个位于重心的十字架来表示，如图6。可以看见，虽然图像中有很多细碎的黑点，但是最后得到的特征块数目并不大（图6中有20个特征块）。这正是因为分合算法具有一定的抗噪声能力。

 

 

 

 

 

第三章  预处理和几何匹配

3.1  预处理

TIFY: inter-ideograph; TEXT-INDENT: 25pt; TEXT-ALIGN: justify">预处理常常有以下几个目的：把图像变成便于后续处理的形式；提高对比度，抑制背景而突出面部特征；去除噪声和其他有害的信息。在我的应用中，把图像变成分合算法能够处理的二值化图像是预处理最重要的功能。其次，希望预处理能够提高二值化前的对比度，把面部器官清晰地突出出来，这样加上合适的二值化方法，就能够保证二值化后的图像中，面部器官以黑色特征块的形式被尽量完整而独立地保留下来。

TIFY: inter-ideograph; TEXT-ALIGN: justify"> 

3.1.1  图像对比度的增强

增强图像对比度有很多种方法，常见的有“S”形变换，直方图均衡化方法等等（参见Kenneth R. Castleman, Digital Image Processing^{（4）TICAL-ALIGN: baseline">）}，都是十分经典而使用很广的，这里就不赘述其具体实现了。简单地说，“S”形变换办法将灰度处于某一范围内（比如灰度很大和很小，或者灰度取值处于中间）的像素的灰度分布差距拉大，从而使对比度提高，这是以牺牲其他灰度范围的对比度为代价的。直方图均衡化则把像素的灰度分布尽量展开在所有可能的灰度取值上，同样能够使对比度提高。

尝试把这些方法应用到人脸定位的预处理中时，效果却不好。原因是这些方法增强对比度没有针对性，结果往往不仅不能突出脸部器官的特征，有时还会造成特征的损失。尤其是直方图均衡化，其结果往往是不能忍受的。

在这样的情况下，我根据实际的需要找到了一种能够突出脸部特征，有针对地增强对比度的方法。

众所周知，边缘提取是特征提取中的一种重要方法。利用边缘提取，可以把图像中灰度发生跳跃变化的地方全部显现出来。脸部器官（尤其是眼睛和嘴）表现为在整个脸的较均匀灰度的背景上的暗团块。所以利用边缘提取可以得到脸部器官的轮廓和内部发生灰度突变的所有边缘。

但是实践证明直接采用边缘提取来提取脸部器官的特征，其效果并不好。边缘提取可类比于微分运算，对于噪声十分敏感，容易造成本不相连的区域粘连起来。

我现在采用的办法是：对输入的原始灰度图像做Sobel运算⁽⁴⁾提取边缘，得到的是灰度的边缘图像（边缘特征用较小的灰度值，即较偏黑色来表示）。将边缘图像按照一定的比例加到原始图像上，再做一个线性变换保证所有灰度值落在0～255的有效灰度范围之内。这样，原始图像中边缘特征就被加强了，脸部器官的特征从而也被加强了，和脸部较为均匀的灰度形成了鲜明的对比。

从图7可以看到原始灰度图像，经Sobel运算提取了边缘之后的灰度边缘图像，以及相加处理完毕的图像（原始图与边缘图比重为5:1）作为比较，使用“S”形变换（参数α^（4）为0.5）和直方图均衡化的结果也显示在图中。

 

还有一种不错的预处理的方法是对灰度分布进行标准化。这是参考了梁路宏, 艾海舟, 何克忠^（7）的方法。预先设定合适的灰度平均值和方差，对输入图像做一个线性变换，使得其灰度平均值和方差都调整得和预先设定值近似相等。这个线性变换的公式是：

  

其中是输入灰度图像，是变换后的灰度图像，x，y是离散化了的像素坐标，W是图像宽度，H是图像高度，是预先设定的标准的灰度平均值和方差，则是输入灰度图像的平均值和方差。可按下式计算：

标准化的效果见图7（标准的灰度平均值为127，方差为10000）。这样做的好处是可以消除部分由于光照条件不同等原因造成的图与图之间的差异，得到一个标准化的输入。而且对于原始输入图像对比度较差的情形，这样做就可以提高对比度。但是实用中发现对于某些情形，这种方法会导致脸部特征的损失。所以目前的算法中没有采用它。

 

3.1.2  二值化

二值化是预处理中最关键的步骤，因为它直接产生能够被分合算法使用的二值化图像。可以说，二值化结果的好坏决定了分合算法结果的好坏，从而决定了整个人脸定位能否成功。

二值化的方法举不胜举，但都可以分为取全局阈值的二值化和取局部阈值的二值化两大类。

取全局阈值的二值化方法中，有一种叫做“组内方差最小化方法”^（8），由Otsu于1979年提出。它的思路是，最好的阈值应该使得被阈值分开的两组的方差的加权和达到最小，其中某组的加权系数就是该组的概率（其实就是该组像素数目占总像素数目的比例）。

设所有像素被阈值t分为两组，灰度≤t的称做组1，灰度>t的称做组2；σ₁²(t)，σ₂²(t)是1组和2组各自的方差；q₁(t)，q₂(t)是1组和2组各自的概率；μ₁(t)，μ₂(t)是1组和2组各自的灰度平均值。定义组内方差σ_W²(t)为这两个组的方差的加权和：

所有像素的平均灰度和方差记做μ和σ²。有下面的关系式成立：

其中σ_B²(t)叫做组间方差，

显然σ²并不随t变化，所以使得组内方差最小的t就是使得组间方差最大的t。寻找这个t的方法是穷举法，即搜索t的每一个可能值，计算出相应的σ_B²(t)，然后找出最大的σ_B²(t)对应的t。

当阈值为t时的各参量已知时，阈值为t+1时的各参量可以用下列递推公式算出，从而避免对每个t都做独立的大量的运算。

，递推起点是q₁(1)=P(1)

，递推起点是μ₁(0)=0

有了q₁(t+1)，μ₁(t+1)和μ₂(t+1)之后，σ_B²(t+1)就可以求出。

这个全局取阈值方法常常可以得到不错的效果。但是对于人脸定位的应用，全局取阈值的方法就不合适了。这是因为人脸各器官相对于面部皮肤的对比度各不相同，比如眼睛就往往比嘴的对比度要高得多。如果对全图像使用同一个阈值，那么无论怎么调整阈值，结果不是眼睛不正确就是嘴不正确，或者两个都不怎么好，捉襟见肘。

我以组内方差最小化方法为基础，构造了一种局域取阈值的方法。其思路是 ：首先对整幅图像用组内方差最小化方法求出阈值，并记录下组间方差。然后把整幅图像划分为m×n个正方形，每个正方形边长都为十几个像素量级。对每个正方形的子图像按组内方差最小化方法求出一个阈值和一个组间方差。最后某一个正方形子图像二值化所使用的阈值是由全局阈值和局域阈值，以全局组间方差和局域组间方差作为权重，再加上事先指定的加大全局阈值权重的因子，综合计算出来。显然全局阈值应该占到绝大部分的比重，因为正方形子图内的像素分部可能与全局分布极不相同，甚至出现全是单一灰度的极端情形。如果仅仅用子图的像素分布来决定子图使用的阈值，那么各个子图取的阈值就会各自相差很大，导致结果图像零乱破碎，无法使用。

全局取阈值二值化方法和局域取阈值二值化方法的结果比较见图8：

 

该局域取阈值二值化方法有边界效应的缺陷。所谓边界效应，是指由于相邻两个正方形子图所取阈值不同，造成边界两边本来灰度变化缓慢的像素有可能一边被二值化成黑色，一边被二值化成白色，形成突变。为了解决这个问题，我曾经构造了另一种局域取阈值二值化方法：对每一个像素，都使用它周围一个边长为十几个像素的正方形区域内的灰度分布来计算一个局域阈值，然后和全局阈值合成这个像素使用的阈值。这样的效果虽然不错，避免了边界效应，但是计算量太大，耗时太长，所以目前的算法中没有采用。这两种局域取阈值二值化方法的结果比较见图8（图中没有明显的边界效应）。