当摄像机从不同视点拍摄同一目标时，所得图像存在一定程度的几何形变，如果目标到摄像机的距离远小于目标大小时，图像的几何形变可用仿射变换近似^[1]。仿射变换下图像的配准有着广泛的应用前景，现有的仿射变换图像配准算法大体可分为局部和全局两类，它们各有优缺点。本文考虑全局的仿射变换配准，这也就是仿射变换参数的恢复问题：给定两幅图像，估计它们之间仿射变换关系中的参数。而全局算法又分为轮廓和区域两类^[2]，轮廓类算法计算量小，但对由多部分构成的目标(如汉字“物”等)无能为力，因此区域类算法受到了更多地关注。

最常用的区域类算法是利用图像矩进行的仿射变换参数恢复。文献^[3]利用交叉权重矩，然而其计算复杂度高；文献^[4]利用高阶矩，其算法对噪声敏感；文献^[5-6]对图像函数进行变换以建立关于仿射变换参数的线性方程组，计算量小，但该类算法对二值图像无能为力；文献^[7]利用图像矩构造非线性方程组，能较方便地解出仿射变换参数，然而该算法所建立的方程组有多组解，需逐个验证，也可能出现无实数解的情形；近来文献^[7]的作者及其合作者们利用与文献^[7]类似的方法做了一系列的工作^[8-9]，但都绕不开非线性方程组多解或无实解的问题。如何设计出计算量小，对灰度、二值图像都适用的参数恢复算法是一个值得研究的问题。

仿射变换的参数恢复问题可通过确定仿射变换前后对应的点来求解。不久前，陈涛等人^[10]提出扩展质心的概念，这种质心与普通质心类似，可保持仿射变换前后的对应关系，结合扩展质心和普通质心可进行图像的规范化、图像识别和仿射变换参数恢复等工作，文献^[10-11]的实验结果也验证了这一点。然而扩展质心利用图像函数的幂来构造，对二值图像无能为力。如何构造对灰度和二值图像都适用、计算量小且能保持仿射变换前后对应关系的类质心点是一个值得研究的问题。

仿射变换的关系式描述为

$\left[ \begin{matrix} {\tilde{x}} \\ {\tilde{y}} \\ \end{matrix} \right]=\left[ \begin{matrix} a & b \\ c & d \\ \end{matrix} \right]\left[ \begin{matrix} x \\ y \\ \end{matrix} \right]+\left[ \begin{matrix} e \\ f \\ \end{matrix} \right]$

(1)

式中,A=$\left[ \begin{matrix} a & b \\ c & d \\ \end{matrix} \right]$是非奇异矩阵，B=$\left[ \begin{matrix} e \\ f \\ \end{matrix} \right]$为平移参数。对矩阵A可进行(奇异值)分解：A=UWV,U、V为正交阵，W为对角阵，因此仿射变换包括旋转、伸缩、平移和斜切变换。

陈涛等人^[10]提出的扩展质心定义为

$\begin{align} & EC{{x}_{\alpha }}=\frac{\iint{{}}x{{I}^{\alpha }}\left( x,y \right)dxdy}{\iint{{}}{{I}^{\alpha }}\left( x,y \right)dxdy} \\ & EC{{y}_{\alpha }}=\frac{\iint{{}}y{{I}^{\alpha }}\left( x,y \right)dxdy}{\iint{{}}{{I}^{\alpha }}\left( x,y \right)dxdy} \\ \end{align}$

(2)

扩展质心E^α(ECx_α,ECy_α)实际上是通过图像函数的幂I^α计算得到的。当α=1时，扩展质心E^α即为普通质心，当α≠1时，E^α为不同于普通质心的扩展质心点，这些扩展质心点与普通质心有类似性质，可保持仿射变换的对应关系，利用扩展质心这一性质可提取仿射不变特征^[10-11]。然而这种扩展质心对二值图像却无能为力(无论α取何值，(ECx_α,ECy_α)都是图像的普通质心)，为此提出一种广义质心，这种质心对二值图像同样适用。

为构造广义质心，将图像原点移至质心，并将直角坐标系转化为极坐标系，极坐标与直角坐标系的变换关系为

$\begin{align} & x=rcos~\theta ,y=rsin~\theta \\ & \tilde{x}=\tilde{r}cos\tilde{\theta },\tilde{y}=\tilde{r}sin\tilde{\theta } \\ \end{align}$

式中，r=$\sqrt{{{x}^{2}}+{{y}^{2}}}$，$\tilde{r}=\sqrt{{{{\tilde{x}}}^{2}}+{{{\tilde{y}}}^{2}}}$，tanθ=$\frac{y}{x}$，$\tan \tilde{\theta }=\frac{{\tilde{y}}}{{\tilde{x}}}$，θ,${\tilde{\theta }}$∈[0,2π)。

定义 对于s,t≥0，令

$\begin{align} & Gx_{s}^{t}=\frac{\int {{\left[ \int {{r}^{s}}f\left( r,\theta \right)dr \right]}^{t}}cos~\theta ~d\theta }{\iint{{}}rf\left( r,\theta \right)drd~\theta }, \\ & Gy_{s}^{t}=\frac{\int {{\left[ \int {{r}^{s}}f\left( r,\theta \right)dr \right]}^{t}}sin~\theta ~d\theta }{\iint{{}}rf\left( r,\theta \right)drd\theta }, \\ \end{align}$

(3)

称点P_s(Gx_s^t,Gy_s^t)为图像的广义质心。

1) 上述广义质心利用变形累次积分计算。图像普通质心利用二重积分定义(如式(2) 中α=1的情形)，通过累次积分计算。而上面的广义质心直接利用变形的累次积分定义(见式(3) )，在对极径方向积分后取幂，然后再对极角积分。

2) 上面定义的广义质心是普通质心的推广。普通质心(式(2) 中α=1的情形)在极坐标系下的形式为

$\begin{align} & EC{{x}_{1}}=\frac{\iint{{}}xf\left( x,y \right)dxdy}{\iint{{}}f\left( x,y \right)dxdy}=\frac{\int [\int {{r}^{2}}f\left( r,\theta \right)dr]cos~\theta d\theta }{\iint{{}}rf\left( r,\theta \right)drd\theta } \\ & EC{{y}_{1}}=\frac{\iint{{}}yf\left( x,y \right)dxdy}{\iint{{}}f\left( x,y \right)dxdy}=\frac{\int [\int {{r}^{2}}f\left( r,\theta \right)dr]sin~\theta d\theta }{\iint{{}}rf\left( r,\theta \right)drd\theta } \\ \end{align}$

也就是说，在式(3) 中s=2,t=1，广义质心对应的就是图像的普通质心，也就是普通质心仅是广义质心的特殊情形。

先将坐标原点移至质心以消除平移，则式(1) 可化为

$\left[ \begin{matrix} {\tilde{x}} \\ {\tilde{y}} \\ \end{matrix} \right]=\left[ \begin{matrix} a & b \\ c & d \\ \end{matrix} \right]\left[ \begin{matrix} x \\ y \\ \end{matrix} \right]$

(4)

此时有如下的定理：

定理 对于s≥0，取t=$\frac{3}{s+1}$，设原图像的广义质心为P_s(Gx_s^t,Gy_s^t)，而经仿射变换后图像的广义质心为${{\tilde P}_s}\left( {G\tilde x_s^t,G\tilde y_s^t} \right)$，则广义质心满足

$\left[ \begin{matrix} G{{{\tilde{x}}}^{t}}_{s} \\ G{{{\tilde{y}}}^{t}}_{s} \\ \end{matrix} \right]=\left[ \begin{matrix} a & b \\ c & d \\ \end{matrix} \right]\left[ \begin{matrix} G{{x}^{t}}_{s} \\ G{{y}^{t}}_{s} \\ \end{matrix} \right].$

(5)

证明 由式(3) 知

$\left\{ \begin{align} & \tilde{r}cos\tilde{\theta }~=arcos~\theta +brsin~\theta \\ & \tilde{r}sin~\tilde{\theta }=crcos~\theta +drsin~\theta \\ \end{align} \right.$

一方面，tan θ=$\frac{y}{x}$，$\tan \tilde{\theta }=\frac{{\tilde{y}}}{{\tilde{x}}}$。因此，$\tan \tilde{\theta }=\frac{ccos~\theta +dsin~\theta }{acos~\theta +bsin~\theta }$，令

$\begin{align} & \alpha \left( \theta \right)=\sqrt{{{(acos~\theta +bsin~\theta )}^{2}}+{{(ccos~\theta +dsin~\theta )}^{2}}}, \\ & \beta \left( \theta \right)=\frac{ccos~\theta +dsin~\theta }{acos~\theta +bsin~\theta } \\ \end{align}$

则

$\tilde{r}=\alpha \left( \theta \right)r,tan~\tilde{\theta }=\beta \left( \theta \right)$

(6)

另一方面

${{(sin\tilde{\theta })}^{2}}=\frac{{{(tan\tilde{\theta })}^{2}}}{1+{{(tan\tilde{\theta })}^{2}}}=\frac{{{(ccos~\theta +dsin~\theta )}^{2}}}{{{\alpha }^{2}}\left( \theta \right)}$

且cos${\tilde{\theta }}$与acos θ+bsin θ,sin ${\tilde{\theta }}$与ccos θ+dsin θ同号，从而

$cos\tilde{\theta }=\frac{acos~\theta +bsin~\theta }{\alpha \left( \theta \right)},sin\tilde{\theta }=\frac{ccos~\theta +dsin~\theta }{\alpha \left( \theta \right)}$

再由式(6) 得${\tilde \theta }$=arctan$\frac{ccos~\theta +dsin~\theta }{acos~\theta +bsin~\theta }$，从而

$\begin{align} & d\tilde{\theta }=\frac{d\left( \beta \left( \theta \right) \right)}{1+{{\left( \frac{ccos~\theta +dsin~\theta }{acos~\theta +bsin~\theta } \right)}^{2}}}= \\ & \frac{adco{{s}^{2}}\theta +adsi{{n}^{2}}\theta -bcsi{{n}^{2}}\theta -bcco{{s}^{2}}\theta }{{{\alpha }^{2}}\left( \theta \right)}d\theta = \\ & \frac{ad-bc}{{{\alpha }^{2}}\left( \theta \right)}d\theta =\frac{\left| A \right|}{{{\alpha }^{2}}\left( \theta \right)}d\theta . \\ \end{align}$

再由${\tilde{r}}$=α(θ)r得到d${\tilde{r}}$=α(θ)dr，因此

$\begin{align} & G{{{\tilde{x}}}^{t}}_{s}=\frac{\int {{[\int {{{\tilde{r}}}^{s}}\tilde{f}\left( \tilde{r},\tilde{\theta } \right)d\tilde{r}]}^{t}}cos(\tilde{\theta })d\tilde{\theta }}{\iint{{}}{{{\tilde{r}}}^{s}}\tilde{f}\left( \tilde{r},\tilde{\theta } \right)d\tilde{r}d\tilde{\theta }}= \\ & \frac{\left| A \right|\int {{[\int {{r}^{s}}{{\alpha }^{s+1}}\left( \theta \right)f\left( r,\theta \right)]}^{t}}(acos~\theta +bsin~\theta ){{\alpha }^{-3}}\left( \theta \right)d\theta }{\left| A \right|\iint{{}}rf\left( r,\theta \right)drd\theta }= \\ & \frac{\int {{\alpha }^{t\left( s+1 \right)-3}}\left( \theta \right){{[\int {{r}^{s}}f\left( r,\theta \right)dr]}^{t}}(acos~\theta +bsin~\theta )d\theta }{\iint{{}}rf\left( r,\theta \right)drd\theta } \\ \end{align}$

令t=$\frac{3}{s+1}$，则有G${\tilde{x}}$_s^t=aGx_s^t+bGy_s^t。

同理有G${\tilde{y}}$_s^t=cGx_s^t+dGy_s^t，命题得证。

该定理说明广义质心与普通质心类似，能够保持仿射变换关系，以一个二值图像的例子说明该定理。图 1(a)(b)分别是仿射变换前后的二值汉字“扬”，图中用圆代表图像的普通质心，三角形代表图像的广义质心(这里s=0.1，从而t=$\frac{30}{11}$)，可以看出，类似于图像普通质心，广义质心可保持仿射变换关系。

仿射变换参数的恢复就是利用两幅满足仿射变换关系的图像f(x,y)与$\tilde{f}(\tilde{x},\tilde{y})$来确定式(1) 中6个参数a,b,c,d,e,f，需要6个独立的条件，只需确定3组满足式(1) 关系的对应点，与扩展质心类似，选择不同参数的广义质心可能在同一直线上，用以下方式确定这3组点：

1) 确定图像的普通质心对P₀(x₀,y₀)和${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{0}},{{{\tilde{y}}}_{0}})$，选择参数s≥0，确定广义质心P_s(x₀,y₀)和${{{\tilde{P}}}_{s}}({{{\tilde{x}}}_{0}},{{{\tilde{y}}}_{0}})$；

2) 以P₀(x₀,y₀)和P_s(x₀,y₀)、${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{0}},{{{\tilde{y}}}_{0}})$和${{{\tilde{P}}}_{s}}({{{\tilde{x}}}_{0}},{{{\tilde{y}}}_{0}})$的连线将图像分割成两部分，确定各部分的质心P₀(x₁,y₁)和P₀(x₂,y₂)、${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{1}},{{{\tilde{y}}}_{1}})$和${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{2}},{{{\tilde{y}}}_{2}})$；

3) 以P₀(x₀,y₀)和P₀(x₁,y₁)、${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{0}},{{{\tilde{y}}}_{0}})$和${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{1}},{{{\tilde{y}}}_{1}})$的连线将图像分割成两部分，确定各部分的质心P₀(x₃,y₃)和P₀(x₄,y₄)、${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{3}},{{{\tilde{y}}}_{3}})$和${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{4}},{{{\tilde{y}}}_{4}})$。

以3组对应点P₀(x₁,y₁)和${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{1}},{{{\tilde{y}}}_{1}})$、P₀(x₂,y₂)和${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{2}},{{{\tilde{y}}}_{2}})$、P₀(x₃,y₃)和${{{\tilde{P}}}_{0}}({{{\tilde{x}}}_{3}},{{{\tilde{y}}}_{3}})$建立方程组确定仿射变换参数。

图 2是图 1中图像经质心与广义质心连线将原图像切片，各切片的质心也标注在图上，由图中可以看出图像经质心与广义质心连线所得切片的质心也保持仿射变换关系，这里只显示一次切片，另一次结果类似。利用这些切片所得质心可将图像恢复，如图 3所示，图 3(a)是利用本文方法恢复的图像，图 3(b)是恢复后图像与原图的差，可看出，恢复的图像与原图像相差不大。

由于计算广义质心时已将坐标原点移至质心，从而平移已消除，因此实验中仅测试式(4) 中矩阵A=$\left[ \begin{matrix} a & b \\ c & d \\ \end{matrix} \right]$的恢复。衡量仿射变换参数恢复效果的公式为

${{e}_{1}}=\frac{\left\| \bar{A}-A \right\|}{\left\| A \right\|},{{e}_{2}}=\frac{\left\| \bar{f}-f \right\|}{\left\| f \right\|}$

(7)

式中，${\bar{A}}$表示恢复的矩阵A，${\bar{f}}$表示恢复的图像f。

以图 1中的二值汉字图像“扬”为例进行仿射变换参数的恢复，图 1(b)是图 1(a)中汉字的仿射变换，仿射变换矩阵为A=$\left[ \begin{matrix} 2.3635 & 1.2860 \\ -0.4168 & 0.3825 \\ \end{matrix} \right]$，选s=0.1,t=$\frac{30}{11}$,利用上述广义质心算法恢复的矩阵为$\bar{A}=\left[ \begin{matrix} 2.3658 & 1.2856 \\ -0.4112 & 0.3857 \\ \end{matrix} \right]$，图 3(a)是利用本文方法恢复的图像，图 3(b)是恢复后图像与原图的差，可以看出，恢复的图像与原图像相差不大，还原后的仿射变换矩阵${\bar{A}}$与正确的A也很相近。

与3.1类似，以图 1中的二值汉字图像“扬”为例，选取相同的仿射变换矩阵A=$\left[ \begin{matrix} 2.3635 & 1.2860 \\ -0.4168 & 0.3825 \\ \end{matrix} \right]$，图 4、图 5分别为文献^[3]，文献^[7]方法对二值汉字“扬”，在相同的变换矩阵下恢复的效果图和恢复图与原图差效果图。

使用文献^[3]所提方法的恢复矩阵为R₁=$\left[ \begin{matrix} 2.2275 & 1.2014 \\ -0.3667 & 0.39717 \\ \end{matrix} \right]$，这与仿射变换矩阵A也很接近，但没有本文算法恢复出的矩阵${\bar{A}}$更接近A，而文献^[7]中方法对应的恢复矩阵为R₂=$\left[ \begin{matrix} 6.1555 & -2.0612 \\ -4.1727 & 3.6806 \\ \end{matrix} \right]$，这与矩阵A相差很大，同时可以从每个算法恢复出的图像与原图像差的效果图看出，本文算法对汉字“扬”的恢复效果要好于文献^[3]和文献^[7]方法。

选用图 6中大小为128×128像素的二值汉字图“扬”作为测试对象，在PC环境( 2.40 GHz CPU，4 GB RAM)，软件Matlab7.0下，对本文算法与文献^[3]和文献^[7]各做一次仿射变换参数恢复耗时比较，结果如表 1 所示,文献^[3]算法的耗时较大，文献^[7] 算法与本文算法耗时相差不大。

为了测试本文算法，与文献^[1]类似以矩阵

$A=l\left[ \begin{matrix} cos~\theta & -sin~\theta \\ sin~\theta & cos~\theta \\ \end{matrix} \right]~\left[ \begin{matrix} a & b \\ 0 & 1/a \\ \end{matrix} \right]$

(8)

来生成仿射变换。仿射变换参数随机选取，选取l∈{0.8,1.2},a∈{1,2},θ∈{0⁰,72⁰,144⁰,216⁰,288⁰，288⁰}，b∈{－1.5,－1,－0.5,0,0.5,1,1.5}因此实验中将每幅图像进行168次变换，并用上面所提算法进行参数恢复。

使用本文算法所得的参数恢复结果与文献^[3]和文献^[7]中的方法所得的参数恢复结果进行比较。表 2，表 3分别给出了实验选取图 6 中的10个二值汉字图像及Coil-20库中20个灰度图像，利用式(8) 对每幅图像各做168种仿射变换，恢复的误差取平均值。由于文献^[7]中利用方程组求解仿射变换参数，会存在方程组无解的情况，本文实验是在剔除无解的情况下对误差取平均值。

由表 2和表 3可以看出，本文算法参数恢复的结果优于文献^[7]，而与文献^[3]的结果相当。文献^[7]利用图像矩构造非线性方程组以求解仿射变换参数，然而该算法所建立的方程组可能有多组解，也可能无实数解，并且所用的矩包括三阶的，然而高阶矩计算不稳定，从而利用该算法不仅有无法实现参数恢复的，即使能恢复的其结果也不佳(文献^[7]所提供的结果都是在大尺寸图像上的实验结果)；文献^[3]基于交叉权重矩进行参数恢复，所用的一阶矩其实就是图像的质心，然而该算法为计算交叉权重矩需要计算图像中每一点的权重，尽管该文算法与本文算法实验结果类似，但其计算量偏大。因此本文算法对二值图像和灰度图像均适用，并且计算量较小。

本文将图像质心推广，提出广义质心的概念，给出了广义质心保持仿射变换前后对应关系的条件，利用广义质心给出了仿射变换参数恢复的算法，该算法对灰度和二值图像均适用，实验结果也表明所提算法优于文献^[3]与文献^[7]所提的算法，但本文考虑全局的仿射变换配准是基于全局的方法，由于利用图像中所有像素点，故计算复杂度比较高，另外基于全局的方法对背景上的噪声点是比较敏感的。由于本文算法是在理想状况下的研究，例如图像背景是比较干净的，但在实际情况中不是这样的，接下来的工作，考虑将该方法运用到实际中的目标。