基于灰度值的匹配方法研究（一）
这在机器视觉领域是一个古老问题，已经被无数的人研究过，提出了各种匹配方案。在OpenCV里也有灰度匹配的函数，只可惜它不支持模板旋转，因此无法应用于实际场景。现在在这里叙述一些我做这个算法后的感悟。供新入门的读者参考。因为一些原因，本系列文章只做流程上的描述和代码片段而不会给出完整的代码。
作为启发性文章，这里说一下灰度匹配是做什么的。灰度匹配是用一个已知的图案作为模板，在运行时图像中查找定位这个模板，一般至少包括模板的位置和角度信息，可以用${ \left( x, y, \theta \right) }$表示。更先进的算法还会输出模板的长宽缩放比例，甚至三维倾斜角等信息。本文不讨论这些更先进的算法。在机械手相关应用中它可以确定物体的位置用于引导机械手抓取东西。而在其它应用中可能只为了确定某个物体的数量或有无。比如用机器视觉检查一个电路板中某个电容有没有正确安装就可以用匹配检测电路板中特定位置有没有出现电容的图案。如下图所示，在一个人像中定位眼睛的位置。图片来自于网络。
本文以基于NCC方法的单目标匹配为例叙述我的实现方法。NCC是一种归一化的相似度衡量方法，非常适用于有相似度阈值的情况。将分为三个部分叙述。第一部分是整体的算法流程；第二部分是相似度计算；第三部分是适应旋转操作的坐标点处理的讨论。
一、算法流程
像其他类型的匹配一样，为了提高效率要对模板和运行时图像做金字塔处理。执行时使用图像金字塔逐层求精得到结果。另外，为了适应模板的旋转，不选择计算模板中的所有像素点，那样速度慢而且旋转模板会导致边缘外出现空白像素，处理起来很麻烦。我的做法是从模板中选择一部分像素点。考虑到随机或等距抽样可能会导致匹配时位置漂移，一个合适的策略是对模板分区，在每个分区中：
- 如果这个分区有边缘则选择边缘最强的点和它附近的两个点；
- 如果这个分区中没有边缘（梯度值普遍较小）则取此分区正中心的点和它附近的两个点。
设置分区的原因是尽可能让选择到的点均匀分布在模板的各个区域，避免点集中在某一个区域里。分区数量可以自定义，我设置的是25个分区。因此总共抽取25×3=75个像素点。下面是我设计的程序的主要流程：
| 参数 |
|
| 生成模板 |
|
| 匹配过程 |
|
| 说明 | 上述参数中的coarseDegree决定了算法的效率。它越大程序运行速度越快。但也不能设置的太大因为还要考虑缩图导致的图像细节的丢失。而fineDegree决定了定位的精度。他越小定位精度越高，但同时算法抵抗形变的能力越差（即一旦运行时图像中的模板有形变容易导致匹配不到）。实际使用中要综合考虑定位精度和抵抗形变的能力。 |
二、相似度计算
NCC的相似度函数公式如下。式中S是运行时图像，T是模板图像，E(*)代表某图像的所有像素平均值，R和C分别是模板的行数和列数，i和j表示模板在运行时图像中的位置：
$${ R(i,j)=\frac{ \sum_{r=1}^{R}\sum_{c=1}^{C}\left| S^{i,j}\left( r,c \right) - E\left( S^{i,j} \right) \right|\overbrace{ \left| T\left( r,c \right) - E\left( T \right) \right| }^{ P\left( r,c \right) } } { \sqrt{ \sum_{r=1}^{R}\sum_{c=1}^{C}\left| S^{i,j}\left( r,c \right) - E\left( S^{i,j} \right) \right|^{2} } \underbrace{ \sqrt{ \sum_{r=1}^{R}\sum_{c=1}^{C}\left| T\left( r,c \right) - E\left( T \right) \right|^{2} } }_{Q} } }$$
上式是计算整个模板所有像素点的公式。其中的P(r,c)和Q只与模板有关，是可以提前计算的。我们的方法只取了模板的一部分像素因此不必所有的r,c都计算。依据该公式典型的计算相似度的函数源代码如下。它是按照人的直觉写出来的，基本思路就是在滑动窗口过程中计算每个位置的相似度，参见下方GrayMatch::execute()函数。
//--------------------------------------------------------------------------------------- // 对NCC模型预处理，提前求出模板自身的相关系数 //--------------------------------------------------------------------------------------- void GrayMatch::preprocessNccParam(TplModel& model) { for (auto& layer : model.samples) /* samples的数据类型是vector<vector<Gray2f>> */ { float avr = std::accumulate(layer.begin(), layer.end(), 0.0f, [](float a, const Gray2f& b) { return a + b.gray; }); avr /= layer.size(); float sum = 0; for (auto& bdr : layer) { bdr.gray = fabs(bdr.gray - avr); sum += bdr.gray * bdr.gray; } sum = sqrtf(sum); layer.push_back({ Point2f(NAN, NAN), sum }); } } //--------------------------------------------------------------------------------------- // NCC计算相似度 //--------------------------------------------------------------------------------------- float GrayMatch::calcScore(const Mat& image, const vector<Gray2i>& bdrs, int x, int y) { float iavr = 0; int count = (int)bdrs.size() - 1; for (int i = 0; i < count; i++) { Point2i pos(bdrs[i].pos.x + x, bdrs[i].pos.y + y); iavr += image.at<uchar>(pos); } iavr /= count; float hcor = 0, lcor = 0; for (int i = 0; i < count; i++) { Point2i pos(bdrs[i].pos.x + x, bdrs[i].pos.y + y); float s = fabs(image.at<uchar>(pos) - iavr); hcor += s * bdrs[i].gray; lcor += s * s; } float ccor = bdrs.back().gray; float sum = hcor / (sqrtf(lcor) * ccor); return sum; } //--------------------------------------------------------------------------------------- // 入口函数大致流程 //--------------------------------------------------------------------------------------- void GrayMatch::execute() { /* 第一层遍历搜索的大致流程 */ /* tplSample是第N层的模板数据，数据类型是vector<Gray2f> */ /* runtime是第N层的运行时图像，数据类型是Mat */ while (每一个旋转角度) { Mat scoreMap = Mat::zeros(runtime.rows, runtime.cols, CV_32FC1); vector<Gray2i> rotatedPoints; rotateGrayBorder(tplSample, angle, rotatedPoints); /* 旋转模板数据 */ ... int miny, maxy; /* 它是滑动窗口y轴的范围 */ int minx, maxx; /* 它是滑动窗口x轴的范围 */ for (int offy = miny; offy <= maxy; offy++) { for (int offx = minx; offx <= maxx; offx++) { float score = calcScore(runtime, rotatedPoints, offx, offy); scoreMap.at<float>(offy, offx) = score; } } ... /* 从scoreMap中找最大相似度的姿势 */ /* 保存全局最优姿势供后续使用 */ } ... /* 根据初始姿势逐层求精 */ }
然而这种方法未必是最高效的写法。从上述代码的GrayMatch::calcScore(...)函数可以看出来，它对图像image中的像素点的访问是跨距的而不是连续访问。因变量pos不是连续变化的。这可能导致较低的缓存命中率，从而降低效率。一般可以通过颠倒循环嵌套顺序来将跨距访问转换成连续访问，大概率可以提高运行效率。代价就是需要分配额外的一个矩阵来暂存中间值，提高一些内存占用。示例如下：
//--------------------------------------------------------------------------------------- // 入口函数大致流程 //--------------------------------------------------------------------------------------- void GrayMatch::execute() { /* 第一层遍历搜索的大致流程 */ /* tplSample是第N层的模板数据，数据类型是vector<Gray2f> */ /* runtime是第N层的运行时图像，数据类型是Mat */ while (每一个旋转角度) { Mat scoreMap = Mat::zeros(runtime.rows, runtime.cols, CV_32FC1); vector<Gray2i> rotatedPoints; rotateGrayBorder(tplSample, angle, rotatedPoints); /* 旋转模板数据 */ ... int miny, maxy; /* 它是滑动窗口y轴的范围 */ int minx, maxx; /* 它是滑动窗口x轴的范围 */ /* 下面就是循环颠倒计算相似度的代码 */ Mat iavr = Mat::zeros(runtime.rows, runtime.cols, CV_32FC1); int count = (int)rotatedPoints.size() - 1; for (int i = 0; i < count; i++) { for (int offy = miny; offy <= maxy; offy++) { const uchar* iptr = runtime.ptr<uchar>(rotatedPoints[i].pos.y + offy); float* optr = iavr.ptr<float>(offy); for (int offx = minx; offx <= maxx; offx++) { optr[offx] += iptr[rotatedPoints[i].pos.x + offx]; } } } Mat hcor = Mat::zeros(runtime.rows, runtime.cols, CV_32FC1); Mat lcor = Mat::zeros(runtime.rows, runtime.cols, CV_32FC1); for (int i = 0; i < count; i++) { for (int offy = miny; offy <= maxy; offy++) { const uchar* iptr = runtime.ptr<uchar>(rotatedPoints[i].pos.y + offy); const float* aptr = iavr.ptr<float>(offy); float* hptr = hcor.ptr<float>(offy); float* lptr = lcor.ptr<float>(offy); for (int offx = minx; offx <= maxx; offx++) { float s = fabs(iptr[rotatedPoints[i].pos.x + offx] - aptr[offx] / count); hptr[offx] += s * rotatedPoints[i].gray; lptr[offx] += s * s; } } } float ccor = rotatedPoints.back().gray; for (int offy = miny; offy <= maxy; offy++) { const float* hptr = hcor.ptr<float>(offy); const float* lptr = lcor.ptr<float>(offy); float* sptr = scoreMap.ptr<float>(offy); for (int offx = minx; offx <= maxx; offx++) { sptr[offx] = hptr[offx] / (sqrtf(lptr[offx]) * ccor); } } ... /* 从scoreMap中找最大相似度的姿势 */ /* 保存全局最优姿势供后续使用 */ } ... /* 根据初始姿势逐层求精 */ }
以上代码经过测试效率提升不是很明显。在第一种直观方法平均耗时约600ms的时候，循环颠倒后平均耗时减少了约20ms。可能的原因是缓冲区（使用的Mat类型变量）数量太多。在其它的较为简单的案例中我实际测试效率提升是较为明显的。尽管我们费了很大力气只取得3.3%的效率提升，但这并非徒劳无功。在将对数据的访问转换为连续访问后这段代码有了向量化的希望，可用SSE/AVX等指令继续优化效率。
三、旋转坐标点中的坐标处理
从模板中选中的坐标点一般是像素的索引，它是整数点${ \left( T_{x},T_{y} \right) }$。还有另外一种处理方法，即把整数坐标加上0.5来变成像素中心点，记为${ \left( T_{x}+0.5,T_{y}+0.5 \right) }$。举个例子如下图所示，O是原点也是模板的旋转中心，其中一个整数像素点是Q(1,0)点，对应的像素中心点是P(1.5,0.5)点。从相机成像机制来说P点更能代表成像时光线作用于光敏传感器的效果。P点旋转180°后四舍五入取整是P’(-2,-1)点，而Q旋转180°后是Q’(-1,0)点。二者是有轻微差别的。从理论上来说用像素中心点的效果会更好一些。
在我的少量的实际测试中，发现用整数点作为模板点坐标进行匹配时，运行时图像中的模板姿势为旋转180°左右时匹配结果确实会有轻微的偏移。而用像素中心点则可以减少匹配结果的误差，不过并不总是这样。由于测试数量不多，所以最终那种处理方式好一点还不确定。有兴趣的可以做大量测试来对比哪种方式更好。
下面是根据上述方法实现的算法效果图。开发环境是VS2015、OpenCV430和Qt5.9，测试电脑配置是AMD A4-9125的CPU，测试用图是130万像素的。图中小窗口模板图里的蓝色的点就是选出来的匹配用的像素点：
经过测试此方法定位精度尚可，且可以抵抗一些诸如模糊、少量噪声或其它类型的形变。上图中仔细看会发现定位稍有偏移，原因是运行时图像有形变。该运行时图像并不是模板中的工件旋转约150°之后的拍的图片，而是模板中工件镜像了之后又旋转拍摄的图片（即从工件的反面拍的图片）。下一篇文章将介绍如何在这个匹配里使用OpenCL加速运算。
