图像对齐算法 - 第二部分

Oleg Krivtsov

4.13/5 (22投票s)

2008年4月21日

CPOL

11分钟阅读

99304

4940

实现并比较了前向组合算法和Hager-Belhumeur算法。

align_demo

引言

图像对齐是一个迭代最小化过程，用于匹配两个图像：模板图像T和另一个图像I。在过去20年中，Lucas和Kanade在1981年提出的图像对齐技术在计算机视觉领域得到了广泛应用。该技术利用图像强度梯度信息来迭代地搜索两个图像之间的最佳匹配。

在《Lucas-Kanade 20 Years On: A Unifying Framework》系列文章中，S. Baker 和 I. Matthews 描述并提供了一个基于Lucas-Kanade技术的图像对齐算法的Matlab实现，包括：

前向加法算法（原始的Lucas-Kanade算法）；
前向组合算法；
逆向组合算法（由Baker、Dellaert和Matthews提出）；
逆向加法算法（由Hager和Belhumeur提出）。

每个图像对齐算法都迭代地最小化模板T和图像I之间的平方差之和，。其中，E(p)是图像相似度函数，p是仿射变换参数向量，x是像素坐标向量，I(x)和T(x)分别是图像I和图像T中像素x的亮度值。

在《图像对齐算法 – 第一部分》中，我们描述了上述算法中的第一个和第三个。在本文中，我们将实现另外两种算法：逆向加法算法和前向组合算法。我们还将修改第一种和第三种算法的代码，使它们能够进行比较，最后，我们将展示哪种算法更好（就速度和收敛性而言），并决定是否可以基于其中任何一种算法构建实时目标跟踪系统。

编译示例代码

要进行代码实验，您需要自行编译。首先，您需要安装Visual Studio。Visual C++ Express Edition也符合我们的要求。

然后，您需要从此处下载并安装OpenCV库。OpenCV是Intel的开源计算机视觉库，我们将它作为我们代码的基础。下载OpenCV_1.0.exe并在您的计算机上运行其安装程序。

最后，您需要按以下方式配置您的Visual Studio。在Visual Studio主窗口中，选择工具 -> 选项...菜单项，然后点击项目和解决方案 -> VC++ 目录树项。

在“显示目录为：”组合框中，选择“库文件”。在下面的列表中添加以下行

C:\Program Files\OpenCV\lib

在“显示目录为：”组合框中，选择“包含文件”。在下面的列表中添加以下行

C:\Program Files\OpenCV\cv\include
C:\Program Files\OpenCV\cxcore\include
C:\Program Files\OpenCV\otherlibs\highgui

就是这样！现在，您可以编译并运行项目了。

项目结构

align_demo项目包含以下文件

文件	它们的内容
auxfunc.h, auxfunc.cpp	包含图像仿射变换函数。
forwadditive.h, forwadditive.cpp	包含Lucas-Kanade（前向加法）算法的实现。
forwcomp.h, forwcomp.cpp	包含前向组合算法的实现。
invadditive.h, invadditive.cpp	包含Hager-Belhumeur（逆向加法）算法的实现。
invcomp.h, invcomp.cpp	包含Baker-Dellaert-Matthews（逆向组合）算法的实现。

它做什么？

当您运行align_demo.exe时，您会看到一个控制台窗口。您需要输入图像变形参数（参数向量p的分量）并按Enter键。

之后，您会看到两个窗口，分别名为“template”和“image”。并且，您会看到一个白色的矩形，表示模板T的初始近似位置。

然后，您会被要求按任意键运行Lucas-Kanade算法。请确保“image”或“template”窗口具有焦点，否则按键将无效。

当Lucas-Kanade算法完成后，您将看到摘要信息，包括计算时间、迭代次数和最终的平均误差。最终的平均误差显示了对齐的质量。“image”窗口中的白色矩形显示了最终的模板位置近似。此外，还会创建一个日志文件。日志文件中的每一行包含迭代次数以及该次迭代的平均误差。您可以在align_demo文件夹中找到日志文件（forwadditive.txt）。我们将使用每个算法创建的日志文件来绘制图表并比较收敛速度。

之后，您将被要求按任意键运行Hager-Belhumeur算法，依此类推。对于每种算法，都会显示摘要信息并创建日志文件。

最后，按任意键退出应用程序。

实现Hager-Belhumeur（逆向加法）算法

我们从Hager-Belhumeur算法开始（您可以在此处找到其原始描述），因为它是一种相当特定的方法。其作者试图实现方法的计算效率，因此他们要求仿射变换W具有特定形式。要求是两个雅可比行列式的乘积可以写成的形式，其中是仅与模板坐标相关的矩阵，而仅与仿射变换参数相关的矩阵。

我们在本文第一部分使用的仿射变换矩阵不具有这种形式。因此，我们需要采用另一种仿射变换。我们采用这种形式：。这种仿射变换可以模拟模板的平面内旋转、平移和缩放。

仿射变换参数向量的形式为，其中是旋转参数，和是平移参数，s 是缩放参数。

一个被仿射变换过的像素的坐标计算方法是仿射变换矩阵W(p)与坐标向量的乘积。

因此，雅可比行列式的计算方法如下：；。最后一个矩阵的逆计算如下：。

经过一些代数运算，这两个雅可比行列式的乘积可以写成要求的形式：。

另一个要求是必须是可逆的。我们雅可比行列式的逆是。

Hager-Belhumeur算法

预计算

计算模板T(x)的梯度。
计算。
计算修改后的最陡下降图像：。
计算修改后的Hessian矩阵：。

迭代

使用W(x;p)对I进行仿射变换，计算I(W(x;p))。
计算误差图像I(W(x;p))-T(x)。
计算。
计算。
计算并更新。

直到。

Hager-Belhumeur算法实现在invadditive.h和invadditive.cpp文件中。

算法的参数包括模板图像T、模板的矩形区域（我们称之为omega）以及图像I。

// Hager-Belhumeur inverse additive method
// @param[in] pImgT   Template image T
// @param[in] omega   Rectangular template region
// @param[in] pImgI   Another image I

void align_image_inverse_additive(IplImage* pImgT, CvRect omega, IplImage* pImgI)
{
    // Some constants for iterative minimization process.
    const float EPS = 1E-5f; // Threshold value for termination criteria.
    const int MAX_ITER = 100;// Maximum iteration count.

实现前向组合算法

为了能够比较前向组合算法和Hager-Belhumeur算法，我们需要在两者中使用相似的仿射变换。让我们采用与Hager-Belhumeur算法相同的仿射变换W 。

前向组合算法

预计算

计算在(x;0)处的雅可比行列式。

迭代

使用W(x;p)对I进行仿射变换，计算I(W(x;p))。
计算误差图像T(x) – I(W(x;p))。
计算图像I(W(x;p))的梯度。
计算最陡下降图像。
计算Hessian矩阵。
计算。
计算参数增量。
更新仿射变换。

直到。

前向组合算法实现在forwcomp.h和forwcomp.cpp文件中。

我们需要在这里存储稍微不同的图像集。我们将存储仿射变换后的图像I(W)、仿射变换后图像的梯度，以及一个包含每个像素处雅可比行列式的8通道图像。

前向组合算法实现在forwcomp.cpp和forwcomp.h文件中。

修改之前实现的算法代码

如果您还记得，我们在本文的第一部分实现了Lucas-Kanade和Baker-Dellaert-Matthews逆向组合算法。现在，我们需要对这些方法的代码进行一些修改，以便与前向组合算法和Hager-Belhumeur算法进行比较。

让我们从Lucas-Kanade算法开始（参见forwadditive.cpp和forwadditive.h）。这里需要修改的地方不多。我们只需将仿射变换矩阵W替换为我们在Hager-Belhumeur算法中使用的那个。并且，由于我们有四个分量的参数向量p，即p=(w_z, t_x, t_y, s)，我们需要使用一个4x4的Hessian矩阵。用于最小化过程的其他矩阵也将具有4x4或4x1的大小。

// forwadditive.cpp; line 50 

H = cvCreateMat(4, 4, CV_32F);
iH = cvCreateMat(4, 4, CV_32F);
b = cvCreateMat(4, 1, CV_32F); 
delta_p = cvCreateMat(4, 1, CV_32F);

我们还将对雅可比行列式和最陡下降图像的计算方式进行一些调整。

// forwadditive.cpp; line 129

// Calculate steepest descent image's element.
float stdesc[4]; // an element of steepest descent image
stdesc[0] = (float)(-v*Ix+u*Iy);
stdesc[1] = (float)Ix;
stdesc[2] = (float)Iy;
stdesc[3] = (float)(u*Ix+v*Iy);

现在，让我们修改逆向组合算法（invcomp.h和invcomp.cpp）。这里，我们不能使用仿射变换矩阵，因为当参数向量p的分量很小时，我们可能会遇到求逆的问题。然而，我们需要使用一个几乎相似的仿射变换来进行比较。让我们将s替换为(1+s)，并使用仿射变换。即使参数向量p的分量等于零，这个矩阵也是可逆的。

这里，我们也有一个4x4的Hessian矩阵，并且需要修改雅可比行列式的计算代码。

比较四种图像对齐算法

由于我们实现的每种算法都使用了几乎相同的仿射变换矩阵和相同的模板，我们可以对它们进行比较。首先，让我们看看哪种算法运行得更快。我们有四个日志文件，其中包含有关我们四种算法收敛性的信息。让我们绘制图表并比较哪种算法收敛得更快。

我们将使用gnuplot工具，它可以从文本文件读取数据并自动绘制图表。要绘制图表，我们将所有日志文件复制到gnuplot的bin目录，然后键入pgnuplot plot.txt，其中plot.txt是一个包含gnuplot命令的脚本（plot.txt也应该复制到gnuplot的bin目录）。

让我们使用Release配置编译align_demo项目，并在我的Pentium IV 3 GHz机器上进行三次简单的实验，使用不同的参数向量p值。您可以在experiments.zip中找到实验结果。

让我们设置p = (0; 5; -5; 1)。这意味着沿X轴平移5个像素，沿Y轴平移(-5)个像素。图表如下所示。X轴表示迭代次数，Y轴表示每次迭代的平均误差值。

让我们设置p = (0.01; 0; 0; 1)。这意味着进行一些平面内旋转。图表放在下面。X轴表示迭代次数，Y轴表示每次迭代的平均误差值。

让我们设置p = (0; 0; 0; 1.05)。这意味着进行一些缩放。图表放在下面。X轴表示迭代次数，Y轴表示每次迭代的平均误差值。

从上面的图表中我们可以看出什么？

在所有三种实验中，算法都能很好地工作，并搜索相似度函数的局部最小值。
在图2中，我们可以看到逆向组合方法会振荡。如果我们添加额外的终止条件，就可以避免振荡。
在图1和图2中，算法收敛很快，并且几乎一半的迭代是不必要的，因为在这些迭代中平均误差值没有减小。所以，我们肯定需要另一个终止条件来避免这些不必要的迭代。

现在，让我们比较一下哪种算法运行得更快。

实验编号	算法名称	计算时间，秒。	迭代次数	最终平均误差
1	前向加法	0.859	100	3.328212
	逆向加法	0.594	100	3.337578
	前向组合	1.078	100	3.367346
	逆向组合	0.469	100	3.341545
2	前向加法	0.875	100	3.397880
	逆向加法	0.579	100	3.395847
	前向组合	1.093	100	3.395613
	逆向组合	0.453	100	7.059829
3	前向加法	0.859	100	5.087051
	逆向加法	0.594	100	24.136780
	前向组合	1.11	100	29.541159
	逆向组合	0.469	100	22.628626

从上面的表格中，我们看到

最快的算法是逆向加法和逆向组合算法。最慢的是前向组合算法。
所有算法都达到了最大迭代次数（100）。也许我们应该增加MAX_ITER常量，即最大可用迭代次数。
在第三次实验中，除了前向加法算法之外，所有算法的最终平均误差都很大。这意味着模板和图像I差异太大，因此需要更多的迭代才能收敛（或者可能根本无法收敛）。如果我们增加MAX_ITER常量，最终的平均误差可能会变小。

结论

在本文中，我们实现了前向组合算法和Hager-Belhumeur逆向加法图像对齐算法。我们对在本文第一部分中已实现的算法进行了少量修改。我们在三个简单的实验中比较了所有已实现算法的性能和收敛速度。

是否有可能基于其中一种算法构建实时应用程序？许多研究人员都表示肯定。例如，J. Xiao等人创建了一个系统，其速度可达约15 FPS。他们的方法与前向组合算法非常相似。

我们尝试将J. Xiao及其同事提出的一些技术应用到我们自己的头部跟踪系统中，您可以在此处找到演示。在我的机器（Pentium IV 3GHz）上，它的运行速度约为15 FPS。使用逆向组合方法有望获得更高的速度。