主动形状模型

机器学习

---

常见使用场景

ASM 模型是一种基于统计形变模型的分割算法。在分割图像时，综合考虑了图像的大小、灰度、大致位置和图像形状等先验知识。它使用从训练样本得到的统计模型作为初始位置，然后通过不断迭代寻找目标图像的最佳轮廓。在迭代过程中通过更新模型的姿态参数和形状参数，使模型逐渐搜索到目标图像的边界。

ASM算法

ASM 是一种基于统计学的模型。通过选取一组训练样本作为训练集，然后用一组标记点来描述训练集中各个样本的轮廓。在算法中，这组标记点将以向量形式使用，称为形状向量。然后对训练集中各样本的形状进行对齐、配准（使得形状消除由于旋转、大小、位置因素造成的差异），使训练集各个样本统一于同一个坐标系。对这些配准后的形状利用主分量分析方法进行降维和统计分析，最后ASM在初始位置的邻近区域进行搜索得到目标形状。由于ASM训练得到的平均形状模型能基本覆盖各种肺部几何外形子空间，具有较高的定位精度，从而可以较准确定位出目标物体。

主动形状模型分为两个步骤：
1. 我们需要训练样本图像得到形状模型和灰度模型。首先，我们需要对训练集样本进行对齐、配准消除因为旋转、拉伸等因素带来的形状差异，然后建立统计形状模型和灰度模型。
2. 我们会使用第一步得到的形状模型和灰度模型去搜索目标图像的轮廓。首先对于测试样本图像，我们需要给出一个初始位置，然后通过灰度模型在初始位置附近进行标记点搜索。当每一个标记点都搜索完毕后会得到新的形状位置，再通过修改主动形状模型的姿态参数进行姿态修正，得到姿态修正的形状向量。最后需要判断形状是否收敛，若否，则返回到标记点搜索这一步进行重新循环直到最后收敛。若收敛，得到的形状就是所求的目标形状。具体流程如下图所示。

建立点分布模型

我们用一组标记点逐一标记训练集中的每一个样本的形状，这样一来，每个样本就可以通过一组标记点来描述，
这组标记点也称为形状向量。Cootes 提出，在ASM 中，标记点可以分为三类：

1. 第一类是解剖学标记点：是由专家指定的，具有某种生物学意义的点，例如心室膜或瓣膜；
2. 第二类是数学标记点：是在目标对象上依靠某些数学或几何学特性确定的点，例如高曲率点，极值点，角点等；
3. 第三类是伪标记点：是在目标对象轮廓上或者标记点直接构造的点，例如前两类点的等距插值点。

下图说明了一般图像中标记点选择方法。

训练样本对齐

ASM 方法通过对样本训练集中的形状进行统计分析，获得形状模型。但是，由于形状在物体的位置，尺寸和角度的影响下，会发生变化。因此，我们需要消除训练集中有关旋转，伸缩和位移因素对形状的影响。我们需要将训练集的每个样本都转化到一个统一的坐标系进行配准，最常用的形状配准方法是Procrustes方法。这个方法通过使训练集中每个样本形状和训练集平均形状之间的距离达到最小来进行配准，其中平均形状可以从整个样本集合中计算得到。假设两个图像的形状向量分别为$x_1$和$x_2$，则两个图像之间的Procrustes 距离为：

$$P^2_d=\sum_{j=1}^n\left[(x_{1j}-x_{2j})^2+(y_{1j}-y_{2j})^2\right]$$
为了使$x_2$对齐$x_1$，我们可通过调整$x_2$ 的旋转角度$\theta$，尺度大小$s$ 和位移$t$ 使其与$x_1$ 间的Procrustes 距离最小，即让(1)式中$E$ 无限接近于0：
$$E=\left(x_1-M(s,\theta)[x_2]-t\right)^TW\left(x_1-M(s,\theta)[x_2]-t\right) \tag1$$
其中
$$M(s,\theta)\begin{bmatrix} x_{2j} \\ y_{2j} \end{bmatrix}=\begin{bmatrix} (s\ cos\theta)x_{2j}-(s\ sin\theta)y_{2j} \\ (s\ cos\theta)x_{2j}+(s\ sin\theta)y_{2j} \end{bmatrix}$$
$W$ 是一个包含每个标记点权重的对角矩阵。权值的选择应该使变化较少的、稳定性高的点权值较大，相反，变化较大的，稳定性差的点使其权值较小，这样做可以更好的保障生成的平均形状的有效性。

点$k$ 权值的计算方法：首先需要计算同一个形状中点$k$ 和点$l$ 之间的距离$R_{ki}$，然后在整个训练集中计算这个距离的方差$V_{R_{kl}}$，然后求点$k$ 上到所有点的距离方
差的和，将这个和求倒数就得到了点$k$ 的权值：

$$w_k=\left(\sum_{l=1}^nV_{R_{kl}}\right)^{-1}$$
通过以上步骤我们已经使训练集中的样本对齐，然后我们需要对对齐后的样本计算平均形状:
$$\overline X=\frac{1}{N}\sum_{n=1}^Nx_i$$
有了以上对两个样本的对齐操作后，就可以对整个训练集采用上述操作，对齐有N个样本训练集的算法如下：
1. 旋转、缩放和平移训练集中每个形状以对齐第一个形状。
2. 计算平均形状。
3. 归一化平均形状的方向、尺度和原点。
4. 重新将每个形状与当前平均形状对齐。
5. 如果过程收敛或者到指定循环次数则退出，否则转到2

形状模型的建立

配准训练集样本后，我们可以很容易求出平均形状，下面对训练集计算协方差$S$:

$$S=\frac{1}{N}\sum_{n=1}^N(x-\overline x)(x-\overline x)^T$$
求出$S$的特征值$p$ 和特征向量$\lambda$ ：
$$Sp=\lambda p$$
如果我们有$p$，其中 $p_i$来自前 $2n$个特征值最大的特征向量。对于任何向量$x$，存在向量$b$ 满足：
$$x=\overline x+Pb$$
这里向量$b$ 的分量表示向量x 在每一个特征向量方向上与平均形状的差异大小。根据PCA我们将选取前$t$个特征值最大的特征向量，来拟合真实的形状，达到降维的目的。

灰度模型的建立

为了使形状模型可以自主地搜索目标图像中肺部的轮廓，我们还需要对训练集中每幅图像标记点附近区域的灰度信息进行采用。首先，对样本中的每个标记点做垂线，再取垂线内外两侧各$k$个像素点的灰度微分，然后对其进行标准化。

为了减少各个样本的整体灰度值差异，我们采样的是标记点的微分而非灰度值。这样我们得到$2k+1$个样本点向量，记作$g_{ij}=[g_{ij1},g_{ij2},...,g_{ij(2k+1)}]$，即第$i$个训练样本形状中第$j$个标记点的轮廓法线向量。如下图所示，图中红线就是第11 个训练样本中第11 个标记点的轮廓法线。

对每个训练样本形状的轮廓法线向量进行标准化，得的第$j$ 个标记点的一个标准轮廓法线向量集。计算关键点的样本平均值$\overline g$ 和样本协方差矩阵$S_g$

$$\begin{align} &\overline g=\frac{1}{N}\sum_{i=1}^Ng_{ij}\\ &S_j=\frac{1}{N}\sum_{i=1}^N(g_{ij}-\overline g_{ij})(g_{ij}-\overline g_{ij})^T \end{align}$$

初始位置的确定

初始位置是由平均形状经过适当的平移，选择，缩放和形变得到。设平移参数为$\vec t=[t_{x1},t_{y1},t_{x2},t_{y2},...,t_{xn},t_{yn}]$，选择和缩放参数分别为$S$ 和$\theta$ ，形变参数为$b$ ，平均形状的形状向量为$\overline X$ ，则初始位置$X_i$ 有：

$$X_i=M(S,\theta)[\overline X+Pb]+t$$
初始化时，设$S=1,\theta=0,b=0,\vec t=\vec 0$，即使用平均形状作为初始位置。

新位置的确定

在前面我们已经得到了模型的初始位置，但是得到的初始位置未必就与目标图像的形状吻合，因此，我们需要调整初始形状的形状参数和姿态参数使其逼近目标形状。在前面我们已经建立了每个标记点的灰度模型，这可以使我们能搜索到每个标记点的最佳位置。

灰度模型是利用从训练集统计出来的标记点在法线方向上的灰度信息来搜索目标样本的的最佳位置，在搜索过程中，标记点在当前位置的法线两侧上各取$m$
个像素点（$m$>$k$）的轮廓法线向量，然后在$2(m-k)+1$个可能位置进行比较，其中，马氏距离最小的那个位置被认为最合适的位置,定义一个函数：

$$f(d)=\left(h(d)-\overline g\right)^TS_g^{-1}\left(h(d)-\overline g\right)$$
这个函数实际上是表示向量$h(d)$和模型均值$\overline g$ 之间的马氏距离，通过最小化$f(d)$就相当于找到一个和训练集的灰度信息最相似的位置。在上面的 $2m+1$ 个像素中分别取其中连续的 $2k+1$ 个像素的灰度值，记为$h(d)$，计算 $f(d)$的值，经过$2(m-k)+1$次计算得到，$2(m-k)+1$个马氏距离，选择其中最小的候选点作为标记点的最佳位置。

计算形状和姿态参数

设当前形状向量为$X$，通过标记点的灰度模型搜索到的新的位置为$X+dX$，我们需要调整当前模型的姿态和形状参数，以便使 $X$ 尽可能的靠近
$X+dX$ ，设姿态参数旋转，缩放，平移分别为$\theta,s,t$，我们可以得到$d\theta,ds,dt$，
其中$X=M(S,\theta)\left[\overline X+Pb\right]+t$ 给出，则$X+dX$ 可以表示为：

$$M(S(1+ds),\theta+d\theta)\left[\overline X+P(b+db)\right]+(t+dt)=X+dX$$

更新形状和姿态参数

由上节我们可以得到姿态参数改变量$d\theta,ds,dt$以及形状参数改变量$db$，我们需要更新这些参数，然后再进行迭代来靠近目标图像的形状。

$$\begin{align} &t=t+dt\\ &\theta=\theta+d\theta\\ &S=S(1+dS)\\ &b=b+db\\ \end{align}$$

算法流程

形状模型和灰度模型建模结束后，便可利用灰度模型收集当前形状下各个标记点的轮廓法线向量，然后选择马氏距离最小的点作为当前标记点的最佳候选点。然后利用形状模型，向最佳候选点构成的形状向量对齐，同时保持形状模型对形状的约束。接着，检查新形状相对于原形状是否收敛。如果没有，则算法结束，搜索完成；否则，继续以新形状的位置寻找最佳候选点，如此反复进行下去，直到定位完成或者算法循环次数达到限定次数。流程图如下图所示：

参考资料

1. http://en.wikipedia.org/wiki/Procrustes_analysis