红外热成像跟踪技术是一种被动式目标检测、跟踪技术，用于对红外视频信号进行目标检测、提取和跟踪。对比度特征鉴别是比较常用的目标提取方法。它无法记忆、识别目标形态特征，在复杂背景下提取效果、跟踪稳定性较差。而模板匹配算法以目标特征数据为模板，在搜索区域里寻找匹配点，即以目标形态特征为判据实现目标检索和跟踪。即便在复杂背景状态下，跟踪灵敏度和稳定度都极高，非常适用于复杂背景下的目标跟踪。
模板匹配算法由于计算量庞大，应用成本较高。经过多方优化、简化后，可用工控机实现实时模板匹配处理。在没有增加成本、耽误工程进度的前提下，增强了复杂背景下的跟踪灵敏度和稳定度，提高了产品的综合竞争能力。为模板匹配算法的低成本应用开辟了一条新路。
本文介绍的模板匹配算法在Ｗｉｎｄｏｗｓ ２０００下用Ｖｉｓｕａｌ Ｃ＋＋编制，可方便地移植到多种操作平台。
１ 模板匹配原理
模板匹配是数字图像处理的重要组成部分之一。把不同传感器或同一传感器在不同时间、不同成像条件下对同一景物获取的两幅或多幅图像在空间上对准，或根据已知模式到另一幅图中寻找相应模式的处理方法就叫做模板匹配。
假设要在搜索区域中寻找与模板图像相关程度最大的位置，可以通过模板匹配来计算两者的相关程度。图１是模板匹配算法的示意图。假设模板（ｂ）叠放在搜索图（ａ）上平移，模板覆盖下的部分记作子图Ｓｉｊ，其中ｉ，ｊ是这块子图的左上角像点在Ｓ图中的坐标。从图１中可得出ｉ，ｊ的取值范围：１≤ｉ≤Ｋ－Ｍ＋１、１≤ｊ≤Ｌ－Ｎ＋１。
衡量模板Ｔ和子图Ｓｉｊ的匹配程度，可用下列两种测度：
Ｄｉｊ＝Ｓｍｎ－Ｔｍｎ（１）
或者
Ｄｉｊ＝（２）
展开前一个式子，有：
Ｄｉｊ＝Ｓｍｎ－２Ｓｍｎ×Ｔｍｎ＋Ｔｍｎ２（３）
|
|
（３）式右边的第三项表示模板的总能量，是一个常数，与（ｉ，ｊ）无关。第一项是模板覆盖下那块子图像的能量，它随（ｉ，ｊ）位置而缓慢改变。第二项是子图像和模板的互相关函数，随（ｉ，ｊ）变化而迅速改变。模板Ｔ和子图Ｓｉｊ匹配时这一项的值最大，因此可以用下列相关函数来反应匹配程度：
Ｒ（ｉｊ＝（４）
或者归一化为：
Ｒ（ｉｊ＝（５）
２ 建立数学模型
２．１ 计算公式
模板匹配算法计算模板和匹配区域的相似程度，以最相似位置为匹配点。由于模板需要在匹配区域上逐次匹配，运算量很大。所以选择匹配公式对整个匹配的效率有极大的影响。
工控机的数据处理能力有限，需要针对红外热成像跟踪技术的特点来简化数学模型，选定计算量最小的计算公式。目标跟踪算法用来确定目标位置，可以用匹配误差的相对大小作为目标判别的依据，误差最小的位置就是目标位置，不需要考虑绝对相似程度。
公式（１）～（５）都能够真实反应模板的相对匹配程度，选择计算量最小、效率最高的公式（１）作为原始数学模型。匹配点位置算法完成整个匹配区域内的最小匹配误差点检索，表示为公式（６）：
Ｄｍｉｎｉｊ＝Ｍｉｎ６
变量Ｋ、Ｌ为匹配区域尺寸；Ｍ、Ｎ为模板尺寸。
２．２ 模板尺寸
模板尺寸对系统性能和计算量的影响不容小觑。模板过大导致动态特性变差过小又会减少目标的特征数据量，降低匹配的敏感程度，增大目标检测难度。实际操作中，模板尺寸设置为３２×１６时效果非常理想。
２．３ 匹配区域
不同的应用环境下，对匹配区域和实时性要求也不尽相同。光电探测设备需要在视频图像采集周期内（２０ｍｓ）完成数据实时处理。由于目标在两场视频图像之间的移动量较小、特征变化不大，匹配区域可以大大缩小。
匹配区域太小会导致目标动态特性变差，过大又会导致计算量大幅度增加，具体选择需要权衡设备参数来决定。由于ＣＣＩＲ制式视频信号是隔行扫描，系统出于实时性考虑，数据以场为单位处理，导致图像比例为２：１状态。为了保持水平、垂直方向的动态特性一致，图像匹配区域也按２：１比例选择。
在满足实时性要求的情况下，选择相对较大的匹配范围，可提高设备的动态特性。从表１实测数据可以看出，选择匹配区域１００×５０点、模板３２×１６点时，动态范围为６９×３５，时间消耗为１３ｍｓ。光电探测设备系统目标动态特性要求处理区域不小于４０×２０点。可见以上选择可以很好地满足动态特性和实时性要求。
３ 数学模型优化方法
数学模型结合选择的模板和搜索区域大小，可以知道模板最佳匹配点计
算公式如下：
Ｄｍｉｎｉｊ＝Ｍｉｎ７
由公式（７）可以看出，程序需要进行大量的循环计算，整体运算量仍然不小，需要进一步优化，减少处理时间。运用如下优化算法进一步减少实际运算量。
３．１ 粗精匹配结合
观察实际模板匹配运算结果可以发现，匹配点附近的匹配误差迅速下降，明显区别于其它位置。针对这一特点，采用粗精匹配结合的算法迅速锁定匹配点大致区域，可大大降低整体匹配次数。
具体实现方法：先跳动着隔几个点进行一次粗匹配，大致框定匹配区域，然后在附近区域逐一检索获得最佳匹配点。运算量可减少到三分之一以下且目标提取效果相当好。
３．２ 限制最大匹配误差
因为只需找到最小匹配误差的位置，不必完整计算每一位置的绝对匹配误差，而以已经计算的最小匹配误差作为最大允许误差。若计算误差大于该最大允许误差，就肯定不是最佳匹配点，可以提前结束计算，进入下一匹配位置的计算；如果匹配完成后仍小于最大允许误差，就用当前误差替换最大允许误差，并把该点作为潜在的匹配位置记录下来。
匹配点和非匹配点的误差常常相差２～３个数量级。经过这种处理后，匹配点后剩余的计算量可以大大降低。
３．３ 乱序匹配
目标出现在匹配区域中的位置不确定。不固定顺序算法可以更快地检索到匹配区域，迅速降低最大匹配误差，减少剩余非匹配点的计算量，降低整体运算量。
针对光电探测设备的实际工作情况，在跟踪状态下，目标位移角速度和角加速度有限，导致目标常处于匹配区域中心附近。选择由中心向周围辐射匹配的方式效果最理想。
４ 程序样本
以下程序样本综合使用了上面的优化算法，成功应用于红外热成像跟踪技术的原理样机，达到了预期效果。
该函数用于图像模板匹配运算，适用于２５６灰度值的黑白图像数据。
Ｄｅａｌ＿ＷｉｔｈＴｅｍｐｌａｔｅＭａｔｃｈｕｎｓｉｇｎｅｄ ｃｈａｒ ｌｐＳｏｕｒｃｅ ＬＯＮＧ ｌＷｉｄｔｈ ＬＯＮＧ ｌＨｅｉｇｈｔｕｎｓｉｇｎｅｄ ｃｈａｒ ｌｐＴｅｍｐｌａｔｅ ＬＯＮＧ
ｌＴｅｍｐｌａｔｅＷｉｄｔｈＬＯＮＧ ｌＴｅｍｐｌａｔｅＨｅｉｇｈｔ
ｕｎｓｉｇｎｅｄ ｃｈａｒ Ｓｏｕｒｃｅ ／／指向待处理图像的指针
ｕｎｓｉｇｎｅｄ ｃｈａｒ Ｔｅｍｐｌａｔｅ ／／指向模板图像的指针
ｉｎｔ ｉｊｍｎ ／／循环变量
ｕｎｓｉｇｎｅｄ ｃｈａｒ ｌＭａｘＷｉｄｔｈ ｌＭａｘＨｅｉｇｈｔ
ｌＭａｘＷｉｄｔｈＥｘａｃｔｌＭａｘＨｅｉｇｈｔＥｘａｃｔ
／／匹配位置
ｕｎｓｉｇｎｅｄ ｌｏｎｇ Ｄ ／／相似误差
ｕｎｓｉｇｎｅｄ ｌｏｎｇ ＭａｘＤ ／／最大允许相似误差
／／粗相关
ＭａｘＤ ＝０ｘ１０００００００ ／／约定最大匹配误差
ｆｏｒ ｊ ＝ ０ｊ ＜ ｌＨｅｉｇｈｔ － ｌＴｅｍｐｌａｔｅＨｅｉｇｈｔ ＋１ ｊ＋＝２
ｆｏｒｉ ＝ ０ｉ ＜ ｌＷｉｄｔｈ － ｌＴｅｍｐｌａｔｅＷｉｄｔｈ ＋ １ｉ＋＝２
Ｄ＝０
Ｓｏｕｒｃｅ＝ｕｎｓｉｇｎｅｄ ｃｈａｒ ｌｐＳｏｕｒｃｅ＋ｌＷｉｄｔｈｊ＋ｉ
Ｔｅｍｐｌａｔｅ＝ｕｎｓｉｇｎｅｄ ｃｈａｒ ｌｐＴｅｍｐｌａｔｅ
ｆｏｒ ｎ＝０ｎ ＜ ｌＴｅｍｐｌａｔｅＨｅｉｇｈｔ ＆＆ Ｄ＜ＭｉｎＤｎ＋＋
ｆｏｒｍ＝０ｍ＜ｌＴｅｍｐｌａｔｅＷｉｄｔｈ ＆＆ Ｄ＜ＭｉｎＤｍ＋＋
Ｄ＋＝Ｓｏｕｒｃｅ＋＋－Ｔｅｍｐｌａｔｅ＋＋Ｓｏｕｒｃｅ＋＋
－Ｔｅｍｐｌａｔｅ＋＋
Ｓｏｕｒｃｅ＋＝ｌＷｉｄｔｈ－ｌＴｅｍｐｌａｔｅＷｉｄｔｈ
ｉｆ Ｄ＜ＭａｘＤ
ＭａｘＤ＝Ｄ
ｌＭａｘＷｉｄｔｈ＝ｉ
ｌＭａｘＨｅｉｇｈｔ＝ｊ
／／精相关
ｌＭａｘＷｉｄｔｈＥｘａｃｔ ＝ ｌＭａｘＷｉｄｔｈ
ｌＭａｘＨｅｉｇｈｔＥｘａｃｔ ＝ ｌＭａｘＨｅｉｇｈｔ
ｆｏｒ ｊ ＝ ｌＭａｘＨｅｉｇｈｔ－２ｊ ＜＝ ｌＭａｘＨｅｉｇｈｔ＋２ ｊ＋＋
ｆｏｒｉ ＝ ｌＭａｘＷｉｄｔｈ－２ｉ ＜＝ ｌＭａｘＷｉｄｔｈ＋２ｉ＋＋
Ｄ＝０
Ｓｏｕｒｃｅ＝ｕｎｓｉｇｎｅｄ ｃｈａｒ ｌｐＳｏｕｒｃｅ＋ｌＷｉｄｔｈｊ＋ｉ
Ｔｅｍｐｌａｔｅ＝ｕｎｓｉｇｎｅｄ ｃｈａｒ ｌｐＴｅｍｐｌａｔｅ
ｆｏｒｎ＝０ｎ＜ｌＴｅｍｐｌａｔｅＨｅｉｇｈｔ ＆＆ Ｄ＜ＭｉｎＤｎ＋＋
ｆｏｒｍ ＝ ０ｍ ＜ ｌＴｅｍｐｌａｔｅＷｉｄｔｈ ＆＆ Ｄ＜ＭｉｎＤｍ＋＋ Ｄ＋＝Ｓｏｕｒｃｅ＋＋－Ｔｅｍｐｌａｔｅ＋＋Ｓｏｕｒｃｅ＋＋
－Ｔｅｍｐｌａｔｅ＋＋
Ｓｏｕｒｃｅ＋＝ｌＷｉｄｔｈ－ｌＴｅｍｐｌａｔｅＷｉｄｔｈ
ｉｆ Ｄ＜ＭａｘＤ
ＭａｘＤ＝Ｄ
ｌＭａｘＷｉｄｔｈＥｘａｃｔ＝ｉ ／／ｘ方向最佳匹配点位置
ｌＭａｘＨｅｉｇｈｔＥｘａｃｔ＝ｊ ／／ｙ方向最佳匹配点位置
/// <summary>
/// 该函数用于对图像进行腐蚀运算。结构元素为水平方向或垂直方向的三个点，
/// 中间点位于原点；或者由用户自己定义3×3的结构元素。
/// </summary>
/// <param name="dgGrayValue">前后景临界值</param>
/// <param name="nMode">腐蚀方式：0表示水平方向，1垂直方向，2自定义结构元素。</param>
/// <param name="structure"> 自定义的3×3结构元素</param>
public void ErosionPic(int dgGrayValue, int nMode, bool[,] structure)
{
int lWidth = bmpobj.Width;
int lHeight = bmpobj.Height;
Bitmap newBmp = new Bitmap(lWidth, lHeight);
int i, j, n, m; //循环变量
Color pixel; //像素颜色值
if (nMode == 0)
{
//使用水平方向的结构元素进行腐蚀
// 由于使用1×3的结构元素，为防止越界，所以不处理最左边和最右边
// 的两列像素
for (j = 0; j < lHeight; j++)
{
for (i = 1; i < lWidth - 1; i++)
{
//目标图像中的当前点先赋成黑色
newBmp.SetPixel(i, j, Color.Black);
//如果源图像中当前点自身或者左右有一个点不是黑色，
//则将目标图像中的当前点赋成白色
if (bmpobj.GetPixel(i - 1, j).R > dgGrayValue ||
bmpobj.GetPixel(i, j).R > dgGrayValue ||
bmpobj.GetPixel(i + 1, j).R > dgGrayValue)
newBmp.SetPixel(i, j, Color.White);
}
}
}
else if (nMode == 1)
{
//使用垂真方向的结构元素进行腐蚀
// 由于使用3×1的结构元素，为防止越界，所以不处理最上边和最下边
// 的两行像素
for (j = 1; j < lHeight - 1; j++)
{
for (i = 0; i < lWidth; i++)
{
//目标图像中的当前点先赋成黑色
newBmp.SetPixel(i, j, Color.Black);
//如果源图像中当前点自身或者左右有一个点不是黑色，
//则将目标图像中的当前点赋成白色
if (bmpobj.GetPixel(i, j - 1).R > dgGrayValue ||
bmpobj.GetPixel(i, j).R > dgGrayValue ||
bmpobj.GetPixel(i, j + 1).R > dgGrayValue)
newBmp.SetPixel(i, j, Color.White);
}
}
}
else
{
if (structure.Length != 9) //检查自定义结构
return;
//使用自定义的结构元素进行腐蚀
// 由于使用3×3的结构元素，为防止越界，所以不处理最左边和最右边
// 的两列像素和最上边和最下边的两列像素
for (j = 1; j < lHeight - 1; j++)
{
for (i = 1; i < lWidth - 1; i++)
{
//目标图像中的当前点先赋成黑色
newBmp.SetPixel(i, j, Color.Black);
//如果原图像中对应结构元素中为黑色的那些点中有一个不是黑色，
//则将目标图像中的当前点赋成白色
for (m = 0; m < 3; m++)
{
for (n = 0; n < 3; n++)
{
if (!structure[m, n])
continue;
if (bmpobj.GetPixel(i + m - 1, j + n - 1).R > dgGrayValue)
{
newBmp.SetPixel(i, j, Color.White);
break;
}
}
}
}
}
}
bmpobj = newBmp;
}
/// <summary>
/// 该函数用于对图像进行细化运算。要求目标图像为灰度图像
/// </summary>
/// <param name="dgGrayValue"></param>
public void ThiningPic(int dgGrayValue)
{
int lWidth = bmpobj.Width;
int lHeight = bmpobj.Height;
// Bitmap newBmp = new Bitmap(lWidth, lHeight);
bool bModified; //脏标记
int i, j, n, m; //循环变量
Color pixel; //像素颜色值
//四个条件
bool bCondition1;
bool bCondition2;
bool bCondition3;
bool bCondition4;
int nCount; //计数器
int[,] neighbour = new int[5, 5]; //5×5相邻区域像素值
bModified = true;
while (bModified)
{
bModified = false;
//由于使用5×5的结构元素，为防止越界，所以不处理外围的几行和几列像素
for (j = 2; j < lHeight - 2; j++)
{
for (i = 2; i < lWidth - 2; i++)
{
bCondition1 = false;
bCondition2 = false;
bCondition3 = false;
bCondition4 = false;
if (bmpobj.GetPixel(i, j).R > dgGrayValue)
{
if (bmpobj.GetPixel(i, j).R < 255)
bmpobj.SetPixel(i, j, Color.White);
continue;
}
//获得当前点相邻的5×5区域内像素值，白色用0代表，黑色用1代表
for (m = 0; m < 5; m++)
{
for (n = 0; n < 5; n++)
{
neighbour[m, n] = bmpobj.GetPixel(i + m - 2, j + n - 2).R < dgGrayValue ? 1 : 0;
}
}
//逐个判断条件。
//判断2<=NZ(P1)<=6
nCount = neighbour[1, 1] + neighbour[1, 2] + neighbour[1, 3]
+ neighbour[2, 1] + neighbour[2, 3] +
+neighbour[3, 1] + neighbour[3, 2] + neighbour[3, 3];
if (nCount >= 2 && nCount <= 6)
{
bCondition1 = true;
}
//判断Z0(P1)=1
nCount = 0;
if (neighbour[1, 2] == 0 && neighbour[1, 1] == 1)
nCount++;
if (neighbour[1, 1] == 0 && neighbour[2, 1] == 1)
nCount++;
if (neighbour[2, 1] == 0 && neighbour[3, 1] == 1)
nCount++;
if (neighbour[3, 1] == 0 && neighbour[3, 2] == 1)
nCount++;
if (neighbour[3, 2] == 0 && neighbour[3, 3] == 1)
nCount++;
if (neighbour[3, 3] == 0 && neighbour[2, 3] == 1)
nCount++;
if (neighbour[2, 3] == 0 && neighbour[1, 3] == 1)
nCount++;
if (neighbour[1, 3] == 0 && neighbour[1, 2] == 1)
nCount++;
if (nCount == 1)
bCondition2 = true;
//判断P2*P4*P8=0 or Z0(p2)!=1
if (neighbour[1, 2] * neighbour[2, 1] * neighbour[2, 3] == 0)
{
bCondition3 = true;
}
else
{
nCount = 0;
if (neighbour[0, 2] == 0 && neighbour[0, 1] == 1)
nCount++;
if (neighbour[0, 1] == 0 && neighbour[1, 1] == 1)
nCount++;
if (neighbour[1, 1] == 0 && neighbour[2, 1] == 1)
nCount++;
if (neighbour[2, 1] == 0 && neighbour[2, 2] == 1)
nCount++;
if (neighbour[2, 2] == 0 && neighbour[2, 3] == 1)
nCount++;
if (neighbour[2, 3] == 0 && neighbour[1, 3] == 1)
nCount++;
if (neighbour[1, 3] == 0 && neighbour[0, 3] == 1)
nCount++;
if (neighbour[0, 3] == 0 && neighbour[0, 2] == 1)
nCount++;
if (nCount != 1)
bCondition3 = true;
}
//判断P2*P4*P6=0 or Z0(p4)!=1
if (neighbour[1, 2] * neighbour[2, 1] * neighbour[3, 2] == 0)
{
bCondition4 = true;
}
else
{
nCount = 0;
if (neighbour[1, 1] == 0 && neighbour[1, 0] == 1)
nCount++;
if (neighbour[1, 0] == 0 && neighbour[2, 0] == 1)
nCount++;
if (neighbour[2, 0] == 0 && neighbour[3, 0] == 1)
nCount++;
if (neighbour[3, 0] == 0 && neighbour[3, 1] == 1)
nCount++;
if (neighbour[3, 1] == 0 && neighbour[3, 2] == 1)
nCount++;
if (neighbour[3, 2] == 0 && neighbour[2, 2] == 1)
nCount++;
if (neighbour[2, 2] == 0 && neighbour[1, 2] == 1)
nCount++;
if (neighbour[1, 2] == 0 && neighbour[1, 1] == 1)
nCount++;
if (nCount != 1)
bCondition4 = true;
}
if (bCondition1 && bCondition2 && bCondition3 && bCondition4)
{
bmpobj.SetPixel(i, j, Color.White);
bModified = true;
}
else
{
bmpobj.SetPixel(i, j, Color.Black);
}
}
}
}
// 复制细化后的图像
// bmpobj = newBmp;
}