前言

一般来说，在计算机视觉任务中，形状匹配可以分为

基于特征（-based）：基于傅立叶变化、特征点、骨架等。 基于亮度（-based）：优先使用像素的灰度值进行形状匹配。

作为一种基于特征的描述符，形状上下文基于两幅图像轮廓采样点之间的相互关系（如距离、梯度）提供了更鲁棒的形状匹配策略。

让我们逐步探索形状上下文算法。

1.轮廓提取（Canny Edge）和轮廓点采样（'s）

轮廓提取

对于轮廓提取，我们使用canny算子，主要步骤如下：

高斯平滑sobel提取梯度幅值和方向非极大值抑制双阈值检测与连接（滞后阈值处理）

这里不再赘述，我们先看一下轮廓采样的步骤。

轮廓采样

对于一般图像，如果采样给定数量的像素，可以采用均匀采样、随机采样、等间隔采样等方法。但对于二维轮廓点集，这些方法就不适用了，因为其分布不均匀。

图1和图2是使用canny算子提取的轮廓。 图3和图4是随机采样100个点得到的采样结果。 可以看出，采样后轮廓信息存在一定损失。 如果我们想要得到如图5、图6这样的采样结果，应该怎么做呢？

一个简单的想法是：我们能否将轮廓点对之间的欧氏距离带入采样步骤，并获得相对均匀（或者说，点对之间充分分散）的轮廓采样结果？ 例如，每次消除最短距离点对中的任意点并消除其与其他点之间的距离关系，重复该过程，直到点对中只剩下所需数量的点集。

这就是主意。 取样流程附如下：

一些技巧：

第一步是随机重新排列轮廓点。 这样可以保证删除的点对的随机性，避免删除时的连锁反应（lab>lbc>lcd>lde。如果每次都删除最小点对，删除的步骤可能会变成e→d→c→b→a ）。 阈值（k=3）主要是针对太多轮廓点预设的多重关系。 如果轮廓点(I)>3×采样点(N)，则直接对轮廓点集进行随机采样，只取前3N个。 可以先对点对距离进行排序，然后将满足条件的右值一一删除（点对中的每个点都没有被删除）（任何左值都可以，只要规则统一）。

下面附上这部分论文样本的源代码：

function [xi,yi,ti]=get_samples_1(x,y,t,nsamp);
% [xi,yi,ti]=get_samples_1(x,y,t,nsamp);
%
% uses Jitendra's sampling method
N=length(x);
k=3;
Nstart=min(k*nsamp,N);
ind0=randperm(N);
ind0=ind0(1:Nstart);
xi=x(ind0);
yi=y(ind0);
ti=t(ind0);
xi=xi(:);
yi=yi(:);
ti=ti(:);
d2=dist2([xi yi],[xi yi]);
d2=d2+diag(Inf*ones(Nstart,1));
s=1;
while s
   % find closest pair
   [a,b]=min(d2);
   [c,d]=min(a);
   I=b(d);
   J=d;
   % remove one of the points
   xi(J)=[];
   yi(J)=[];
   ti(J)=[];
   d2(:,J)=[];
   d2(J,:)=[];
   if size(d2,1)==nsamp
      s=0;
   end
end

2、形状上下文（Shape）直方图矩阵的构建及其相似度度量方法

对于轮廓中的点对，哪些信息有助于形状匹配？

距离、角度等…

那么对于每个轮廓点，我们可以构造一个以其为中心的对数极坐标系（包括12个角度区域和5个距离区域，如下图c所示），然后将其周围的轮廓点映射到每个区域内，然后统计落在每个区域的轮廓点的数量，最后进行归一化（即除以落在所有区域的轮廓点的数量）。 这样我们就生成了一个形状上下文直方图矩阵。 如下图所示，我们可以直观地看到相似度de>相似度df。

一些技巧：

使用极对数坐标的好处主要是使描述子对相邻采样点更加敏感，强化局部特征。 首先获得点对的距离矩阵，然后除以距离均值，以使形状上下文描述符对缩放不敏感。 然后将距离矩阵映射到对数域并与阈值进行比较。 由于对数函数是单调变换，因此对数域中的五个距离阈值可以映射为欧氏距离（0.125、0.25、0.5、1.0、2.0）。这里大于2.0的点被认为距离目标点太远，因此不考虑）形状，然后将距离矩阵与阈值进行比较，以去除对数操作的时间成本。 如下所示：

% use a log. scale for binning the distances
r_bin_edges=logspace(log10(r_inner),log10(r_outer),5);

计算角度矩阵的方法也如上图所示，比较简单。

解出距离矩阵和角度矩阵后，就可以进行区域映射并统计每个区域的点数。 如果我们在这里采样100个轮廓点，就会生成一个大小为100×60的矩阵。

不要忘记进行标准化。

至此，形状上下文直方图矩阵就构造完成了。 下面我们使用卡方统计量 χ2 来计算点输入图和点模板图之间的相似度。

经过前面的计算，我们得到了输入图像和模板图像的N(100)×K(60)形状上下文直方图矩阵。 对于输入图像和模板图像中的每个点，都有一个K维直方图向量。 ，分别记为g(k)和h(k)，代入卡方公式计算相似性测度矩阵：

此时我们得到了关于形状上下文的相似性度量矩阵（N×N）。

这里解释一下为什么是N×N矩阵：输入图像和模板图像都采样了N个轮廓点（也可以采样不同数量的轮廓点，然后通过我们后面会讲的伪点策略进行匹配） ），然后对于两个图片中的每一个点对，将g(k)和h(k)带入卡方公式，就可以计算出一个Cs。 这样总共生成了N×N个矩阵。

3、局部外观（Local）描述矩阵的构建及其相似度度量方法

根据形状上下文直方图矩阵，我们获得全局信息，即全局点到目标点的方向矩阵，类似于不同方向和角度的轮廓上点集的密度。 接下来，我们考虑添加另一个局部外观信息来约束局部特征。

有很多本地外观信息，比如

局部方向（局部）局部颜色直方图（局部颜色）局部纹理（局部）局部拓扑（局部）等等…

我们在这里考虑使用局部方向（local）作为形状上下文的局部外观信息。

对于输入图像和模板图像，我们计算对应灰度图像的梯度Gx和Gy，然后找到轮廓点对应的梯度值并找到切线角θ1、θ2。 这里我们引入切向角相异（角度）函数作为测量局部方向的基础：

为什么可以这样定义呢？ 我们看下面的图片。 假设有一个半径为1的圆，那么对于θ1和θ2，通过余弦和正弦将它们变换到欧氏空间，然后计算它们的欧氏距离，就可以衡量两个角度的相似度了。

一些技巧：

对切​​向角相异函数进行变换，可以得到简单的计算公式：

最后一步是对两个相似性测度矩阵进行加权求和并输出总的相似性测度矩阵：

Cs对应形状上下文（Shape），CA对应局部外观（local），论文中给出β=0.1。

4.匈牙利()策略

匈牙利匹配是一种常见的二分图匹配，复杂度为O(N3)。 常用的场景包括：

求最大匹配问题就是求增广路径的过程：

给定一个匹配关系，找到使匹配数量最大化的匹配。 如下图所示，最大匹配为3个，即L1-R2、L2-R3、L3-R1。 我们不会详细介绍该方法的实现（深度或广度优先匈牙利）。 如果你有兴趣，可以尝试一下。

寻找最短路径的问题：给定N对N个点之间的距离，找到使总距离最小的点对匹配关系。

这里我们将之前的相似度度量矩阵抽象为距离矩阵。 使用匈牙利算法的目的是找到使总损失最小化的点对匹配。 具体步骤如下：

变换距离矩阵Aij→Bij，使得Bij的所有行和列都有0次尝试分配。 如果得到最优解，则返回损失最小。 否则，进入(3)制作覆盖所有0个元素的最小直线，并转至(4)增加距离矩阵。 多余0个元素，转(2)

下面给出两个例子：

如下图所示，给出了4-4个点对的距离矩阵。

对应的匹配关系如下：

如下图所示，给出另外4-4个点对的距离矩阵。

让我们思考一下为什么我们可以这样做：

由于距离矩阵的这种对应关系，无论我们逐行（列）减去对应行（列）的最小值，都不会影响距离排序。 参考例2，在加0的步骤中，我们取出未被直线覆盖的矩阵，并从每个元素中减去最小的元素e。 当我们把这个矩阵放回去时，为了保持距离排序关系，我们需要将位置14和44对应的0元素减去e，或者将e添加到位置24和34。

下面的链接是匈牙利算法在线测试的链接：

经过匈牙利算法，找到两张图片（digit-0，digit-3）的对应关系。 我们用白线连接对应关系：

从上图我们可以直观地发现大部分的点对关系都已经找到了。 但仍然有一些（中间）点突然匹配。 这是我们真正不想要的。 我们是否可以添加一个阈值，使得这些匹配损失非常大的点对不被匹配？

下面我们介绍一下伪点匹配。

5. 伪点匹配

通过添加伪点，

解决两幅图像轮廓点不同的情况。 消除一些匹配损失过多的点对匹配（异常值检测 – ）

我们来仔细看看如何实现伪点匹配。

如果距离矩阵不是N×M（M>N），如图1所示：

异常值检测和去除。 下左图显示了4-4个点对的距离矩阵。 如果不添加伪点，得到的最小距离如右图。

接下来考虑添加伪点的情况。

这意味着对于A点，如果此时没有损失小于阈值的匹配点，则会分配一个伪点。

下图展示了伪点匹配对匹配结果的影响。 左边是没有添加伪点的情况，右边是添加伪点的情况：伪点数量为轮廓点的0.25，损失为0.25。

添加伪点使得轮廓对噪声更加鲁棒，但这也意味着对于下一次迭代（经过TPS变形的图像，TPS将在下一节中介绍），形状上下文矩阵和相似度矩阵将受到外观的影响异常值。 并做出相应的改变：即不再考虑异常值的贡献。 同时，我们仍然需要获取这些离群点的形状上下文，因为这些离群点在下一次变形后可能会变成内部点（ ）。 详细内容请参见源代码。 6． 用于轮廓点拟合的薄板样条插值（Thin Plate）

通过匈牙利匹配和伪点匹配策略，我们从2N个轮廓点中获得了S个匹配度较高的匹配点对。 我们提出了另一个想法：我们是否可以使用这些S匹配点对来使输入轮廓变形，使其更接近模板轮廓？ 那么对于变形轮廓和模板轮廓，如何再次找到匹配度高的匹配点对呢？

因此，我们引入了薄板样条插值（TPS）函数来使输入轮廓变形，以最小化总弯曲能量。

薄板样条插值是一种常见的插值方法，常用于二维图像配准。 通过将N个主要点（Point）变形到对应的位置，得到对应的位移场，然后计算其他区域点的变换后的坐标。 如下所示：

对于一般插值方法，插值的主要约束是：

给定一个插值点(rj,tj)，尝试构造一个连续平滑函数y=f(x)，使得函数曲线(曲面)经过所有插值点，即f(rj)=tj。 因此目标函数可以写成以下形式：

对于多谐波样条插值 ( )，添加基于梯度的平滑项（正则化项）：

薄板样条插值是多调和样条插值的特例 (m=2)

我们将正则化项中的λ称为正则化参数，其积分部分称为弯曲能，即：

我们可以将薄板样条插值函数想象成一块弯曲的薄钢板，将其穿过顶点并找到使弯曲能量最小的插值函数。

其中 Φ( r ) = r2lnr 为 TPS 的径向基核函数，其二维图如下：

文献[1]证明，上述f(x)函数是使弯曲能最小化的插值函数。 可以看作是非变形项/线性项(a0+a1x+a2x)加上变形项(后一部分)。 线性项定义了所有控制点的平面的最佳匹配，并且可以被视为最小二乘拟合，而线性项可以被视为与控制点提供的曲率率相关。 如下图所示，当正则化参数λ越大时，越接近最小二乘拟合。

由于插值函数有 N+3 个条件，并且我们有 N 个（点对）约束，因此我们再添加 3 个约束，如下所示：

这三个条件与插值条件一起可以写成以下矩阵形式：

然后将正则化参数 λ 带入解中（在文章最后的公式推导中），我们可以得到

终于得到了

弯曲能量 Dbe =mean(diag(wT(A)w)) 变形损失 (Cost) = log(svd(s(1)/s(2)), 其中 s = svd(a). 后轮廓点的坐标转型。

下面附上论文样本中TPS/的源代码以及注释：

function [cx,cy,E,L]=bookstein(X,Y,beta_k);
% [cx,cy,E,L]=bookstein(X,Y,beta_k);
%
% Bookstein PAMI89
N=size(X,1); 	% 轮廓点集1
Nb=size(Y,1); 	% 轮廓点集2
if N~=Nb 	% 判断轮廓点集是否相同
   error('number of landmarks must be equal')
end
% compute distances between left points
r2=dist2(X,X) 	% 求距离
K=r2.*log(r2+eye(N,N)); % add identity matrix to make K zero on the diagonal 求A矩阵（径向基部分）
P=[ones(N,1) X]; 	% 构造B矩阵
L=[K  P
   P' zeros(3,3)];	% 构造L矩阵 (Lc = V)
V=[Y' zeros(2,3)]; 	% 构造V矩阵
if nargin>2
   % regularization
   L(1:N,1:N)=L(1:N,1:N)+beta_k*eye(N,N); %如果有正则化项
end
invL=inv(L); 		% 求L的逆
c=invL*V';		% 求c矩阵 这里c = [c a]T
cx=c(:,1); 		% x方向
cy=c(:,2); 		% y方向
if nargout>2
   % compute bending energy (w/o regularization)
   Q=c(1:N,:)'*K*c(1:N,:);
   E=mean(diag(Q)); 	% 求解弯曲能量
end

在多次迭代过程中，获得TPS变形图像和模板图像之间的形状上下文，然后使用匈牙利匹配来查找点对关系。 如下图（左边是变形过程以及各个loss的输出值，右边是匈牙利算法得到的匹配点对）：

7.损失函数计算和k-NN分类器进行分类

经过多次迭代（n=5）后，我们可以计算总损失来确定输入图像与模板图像的相似度：

形状上下文距离 Dsc：

形状上下文距离是通过对称求和得到的，即对于总相似度度量矩阵，加权的每行最小损失值除以列数和每列最小损失值除以行数。 公式如下：

其中 T 是 Q（输入图像）的 TPS 变换。

外观距离Dac：

通过高斯窗函数，计算TPS变换后P（模板图像）中轮廓点与Q（输入图像）中轮廓点的亮度差平方平均值作为中心和外围。

弯曲能量距离Dbe：

弯曲能量是通过薄板样条插值函数得到的值。

最后，使用经验公式来查找模板和输入图像之间的相似度。 公式如下：

这种相似度也可以作为数字识别的k-NN分类器的分类指标。

八、总结

最后给出形状上下文（Shape）的过程：

并且，形状上下文论文的地址和实现

如果想更深入的了解薄板样条插值及推导过程，这里有薄板样条插值源码论文（1篇）和相应公式的推导（2篇）的下载链接

9.一些改进

[1] Kent, JT 和 KV (1994a)。 链接和薄板。 见： 、 和 ： a to Peter（FP Kelly 编辑），第 325-339 页。 约翰威利父子有限公司。 第 282 页、287 页、311 页

[2] S、森 G、马利克 J。 形状 [C]// IEEE 基于图像和视频。 IEEE，2000。

[3] S，马利克，J，J。 [C]// IEEE 上的 . IEEE，2001。

[4] ，S，Malik，J，，J.形状与使用形状[J]。 IEEE on & ，2002，24(4)：509-522。