计算机视觉之统计形状模型——Active Shape Model

ASM（Active Shape Model）由Cootes于1995年在论文[1]中提出。

ASM是对图像中的shape进行建模的可变模型（deformable model），得到的形状模型既可以用来分析新的形状（拟合模型到新形状，见第4节），也可以用于生成形状（在给定图像中搜索形状，见第5节）。

首先需要给出shape的定义。给定一幅图片，对感兴趣物体标注landmark points，以人脸为例，就是选择人脸上的关键点，如下图所示：



得到点集合{(x1,y1),(x2,y2),⋯,(xn,yn)}，那么这些点在图片中对应的shape就可以表示成2n维的向量x⃗ ：

x⃗ =(x1,y1,x2,y2,⋯,xn,yn)T

可以看出每个shape都是R2n空间中的一个点。并且我们规定对一个形状进行相似变换（Similarity transformation），即旋转R、缩放S、平移T（RST），并不会改变该形状。当然也可以定义形状在其他变换操作下是不变的。

假设现在有N个已经标注了landmark点的训练数据，也就是说有了N个形状{x⃗ 1,x⃗ 2,⋯,x⃗ N}。

1、对齐

为了能比较不同shapes，需要将shapes对齐。这是因为不同的形状大小角度等都可能不同，比如一个人脸特写与一张全身照中的人脸，显然，这两个人脸对应的形状是没有可比性的。

ASM中的对齐算法是基于Procrustes Analysis的。Procrustes Analysis在对齐形状时以最小化对齐后的形状和平均形状之间的距离之和为目的：

D=∑i|x⃗ ′i−x⃗ ¯|2.(1)

这个优化目标有意思，首先式中的x⃗ ′i是对齐后的第i个形状，x⃗ ¯是对齐后的形状的均值。第一项是未知的，而第二项又取决于第一项。假设x⃗ ′i=T(Θi;x⃗ i)，这个函数姑且称为对齐函数吧。那么每个形状都对应一个参数Θi，上述优化问题就是求使D最小的Θ={Θi,i=1,⋯,N}：

Θ∗=argminΘD=argminΘ∑i|T(Θi;x⃗ i)−1N∑jT(Θj;x⃗ j)|2(2)

这个优化问题有解析解，但是一般采用下面的迭代算法进行对齐：

平移每个形状使其重心位于原点

选择第一个形状x⃗ 1，对齐进行缩放x⃗ 1以使|x⃗ |=1，以缩放后的形状作为初始平均形状x⃗ 0¯

对齐所有形状到平均形状

计算对齐后的所有形状的平均形状x⃗ ¯

对x⃗ ¯进行约束：让x⃗ ¯对齐x⃗ 0¯并且缩放到|x⃗ ¯|=1

如果不收敛，返回到(3)

说明：步骤(5)对平均形状进行约束，是为了确保有唯一解（如果不进行约束的话，平均形状可能变得非常离谱）。另外判断是否收敛是看平均形状的变化是否超过某一幅度。

1.1、对齐操作

步骤(3)中的对齐操作可以有多种方式。比如可以移动每个形状使其以原点为中心，缩放使其模为1，旋转使得D最小，这种方式下D的参数只有旋转角度。

文献[1]中介绍的方法是对每个形状进行相似变换Ts,θ,tx,ty(x⃗ )以对齐均值。假设两个形状x⃗ 1和x⃗ 2，这两个形状以原点为中心，现在x⃗ 1处理做相似变换使其对齐x⃗ 2，相似变换的参数通过最小化

E=|Ts,θ,tx,ty(x⃗ 1)−x⃗ 2|2(3)

来确定。相似变换的表达式为

Ts,θ,tx,ty(x⃗ )=s[cosθsinθ−sinθcosθ]x⃗ [⋅]+[txty]

变换得到：

Ta,b,tx,ty(x⃗ )=[ab−ba]x⃗ [⋅]+[txty](4)

其中x⃗ [⋅]表示遍历x⃗ 中的每一个点。

现在来求解这个优化问题。令E对每个参数的偏导数为0，得到下面的公式：

a∑i=12n[(xi1)2+(yi1)2]+tx∑i=12nxi1+ty∑i=12nyi1=b∑i=12n[(xi1)2+(yi1)2]+ty∑i=12nxi1−tx∑i=12nyi1=a∑i=12nxi1−b∑i=12nyi1+tx=a∑i=12nyi1+b∑i=12nxi1+ty=∑i=12nxi2xi1+∑i=12nyi1yi2∑i=12nxi1yi2−∑i=12nyi1xi2∑i=12nxi2∑i=12nyi2(5)

假设x⃗ 1已经平移到以原点为中心，那么∑2ni=1xi1=∑2ni=1yi1=0，于是有

txtyab=∑i=12nxi2=∑i=12nyi2=∑2ni=1xi2xi1+∑2ni=1yi1yi2∑2ni=1[(xi1)2+(yi1)2]=x⃗ 1⋅x⃗ 2||x⃗ 1||2=∑2ni=1xi1yi2−∑2ni=1yi1xi2∑2ni=1[(xi1)2+(yi1)2]=∑2ni=1xi1yi2−∑2ni=1yi1xi2||x⃗ 1||2(6)

如果x⃗ 1没有移动到以原点为中心，那么需要求解复杂的方程组(5)。

最终得到：

s=a2+b2−−−−−−√θ=tan−1ba

注意文献[1]在式(3)中加了一个权重矩阵，以区别不同的点对距离的影响程度，感兴趣的可以阅读[1]。

2、对形状的变化进行建模

对齐后的x⃗ 在2n维空间中构成一个cloud，这个cloud对应这一个分布。如果能得到这个分布，那么我们不仅可以生成和训练数据相似的新形状，还可以分析给定的一个形状是否和训练数据一致。

2.1、PCA

首先，我们采用PCA降维技术对数据进行降维以便处理。PCA降维过程如下：

首先计算数据均值：x¯=1N∑i=1Nx⃗ i

然后计算数据的协方差矩阵：S=1N−1∑i=1N(x⃗ i−x¯)(x⃗ i−x¯)T

计算S的特征值λi和对应的特征向量p⃗ i，并且按特征值递减顺序进行排列，即λi≥λi+1

用S的前t个特征值对应的特征向量构成矩阵Φ（2n×t），那么数据x⃗ 可以近似表示成x⃗ =x¯+Φb⃗ (7)

2.2、t的确定

我们知道，每个特征值λi给出了训练数据在相应特征向量方向上的方差，那么训练数据的总方差为V=∑λi.t可以通过下式来确定：

∑i=1tλi≥fV

其中f表示希望保留多少可变性。

2.3、分析b⃗

由式(7)可知，不同的b⃗ 对应不同的形状x⃗ 。我们认为和训练数据一致的形状是合理的，那么自然b⃗ 也存在一个合理范围，就是说b⃗ 的取值不是任意的，需要约束到合理范围内。我们假设b⃗ 服从分布p(b⃗ )，p(b⃗ )可以根据训练数据进行估计。我们认为当p(b⃗ )≥pt(8)时，形状是合理的。

假设bi是独立的并且服从高斯分布，那么有logp(b⃗ )=−0.5∑i=1tb2iλi+const.(9)

我们可以对bi进行硬约束−3λ√i≤bi≤3λ√i.(10)

也可以根据公式(8)并且结合公式(9)得到另一种约束：

(∑i=1tb2iλi)≤Mt(11)

其中左边是b⃗ 和均值之间的Mahalanobis距离，右边的Mt可以通过卡方分布得到。

假设b⃗ 服从高斯分布可以解决大部分的问题，但是无法表示非线性的形状变换，比如只对形状的一部分进行旋转得到的训练数据。更具体地，我们通过旋转在一个正方形中的三角形得到一组训练数据，生成的样本如下图所示。



对训练数据应用PCA技术，发现有两个主要成分。我们把数据映射到二维空间中得到b⃗ 的分布如下图所示：



显然，b⃗ 不是高斯分布。

这时可以用混合高斯模型GMM对p(b⃗ )进行建模：

p(b⃗ )=∑j=1mwjG(b⃗ ;μj,Sj)

3、为当前模型形状x⃗ ′找到最近的可能模型形状x⃗

这个标题有点拗口，现在来解释一下为什么要给形状加定语“模型”。对于给定的一组训练形状数据{X}，把每个形状看作2n维空间中的一个点，可以想象得到这些点在2n维空间中分布比较散乱，这就是因为形状还没有对齐。假设经过对齐（相似变换）后得到的形状集合为{x⃗ }，再在2n维空间中找到对应的点，那么这些点就比较紧凑了。我们为了区分这两个形状，把X称为图像形状，把x⃗ 称为模型形状。这里处理的是模型形状。

而根据在PCA小节的介绍，每个模型形状可以由公式(7)表示，于是有

x⃗ ′=x¯+Φb⃗ ′

可以得到：

b⃗ ′=ΦT(x⃗ ′−x¯)

我们在2.3节中提到过当p(b⃗ )≥pt时，认为形状是合理的。而当p(b⃗ )<pt时，形状是不合理的，但是我们希望找到离形状最近的合理形状。这时无论b⃗ 是高斯分布还是混合高斯分布，我们对p(b⃗ )求梯度，然后沿着上坡方向移动直到超过阈值pt。

4、用模型拟合新的形状(Fitting a model to new shape)

现在我们讨论怎么求得最优pose参数{s,θ,Xt,Yt}和shape参数b⃗ 以匹配给定的图像形状Y。

在前面的的形状对齐过程中，我们是对图像形状X做相似变换T，映射到模型形状x⃗ 。那么也可以用一个相似变换T′把模型形状x⃗ 映射为图像形状X，并且T′可以由T得到

T′=T−1

T−1=1s[cosθ−sinθsinθcosθ]([xy]−[txty])

针对这一节的问题，我们通过优化问题

min||Y−Ts,θ,tx,ty(x¯+Φb⃗ )||2

来求得最优参数。

下面给出优化问题的迭代法：

初始化b⃗ =0

得到模型形状：x⃗ =x¯+Φb⃗

寻找最优参数{s,θ,tx,ty}使得最好地把x⃗ 匹配到Y：即min||Y−T(x⃗ )||2，这一步骤参考1.1节

通过逆变换把Y映射为模型形状：y⃗ =T−1(Y)

更新参数：b⃗ =ΦT(y⃗ −x¯)

对b⃗ 进行约束

判断是否收敛，如果不收敛则返回步骤2

判断是否收敛是通过检查参数是否变化。

5、图像搜索——ASM

给定一张图片怎么找到其中的形状呢？

初始化x⃗ =x¯， b⃗ =0，并在图像中初始化一个形状X

检查X每个点附近的区域，找到最好的新的形状X′

更新参数s,θ,tx,ty,b⃗ 以拟合新形状X′：第4节

返回2直至收敛

这里需要着重介绍怎么根据当前形状X找到更好的形状X′。

对于以边缘作为形状的情况，可以沿着模型边界的法线方向移动，直至到达图像上梯度较高的点。如下图所示：



但是有些情况下形状不一定就是边缘，可能表示其他的结构。这时，可以通过对形状附近的图像结构进行建模，在搜索时，找最匹配所建模型的点。

常见的是对灰度结构进行建模。对训练形状Xi中的一个模型点沿法线方向在两侧各取k个点，于是每个模型点都得到2k+1个点，用灰度值构成向量g⃗ i。假设有N个训练数据，对一个模型点我们可以得到N个向量{g⃗ i}。我们假设g⃗ 服从高斯分布，并且可以得到它的均值g¯和协方差Sg。这样对一个模型点我们得到一个统计模型。

给出一个新的样本g⃗ s，Mahalanobis距离：

f(g⃗ s)=(g⃗ s−g¯)S−1g(g⃗ s−g¯)T

和g⃗ s来自于这个分布的log概率呈负相关性。

在搜索时，我们检查当前形状的模型点沿着模型边界的法线方向两侧各m(m>k)个点，这样可以得到2(m−k)+1个g⃗ （自己在纸上画一画就知道为什么是这么多），找到使得Mahalanobis距离最小的g⃗ 作为新的模型点。

6、实现

Stasm：作者Stephen Milborrow，他研究生和博士生阶段一直在研究ASM，于2016年博士毕业。在另一篇博客中简单介绍了Stasm源码。

参考：

[1] Cootes, T.F., Taylor, C.J., Cooper, D.H., Graham, J.: Active shape models—their training and application. CVIU 61, 38-59 (1995)

[2] Cootes, T.F., Taylor, C.J.: Technical Report: Statistical Models of Appearance for Computer Vision. The University of Manchester School of Medicine (2004), www.isbe.man.ac.uk/∼bim/refs.html

[3] Youtube视频：https://www.youtube.com/watch?v=53kx_czs7Es