【工程应用一】 多目标多角度的快速模板匹配算法（基于NCC，效果无限接近Halcon中........)

       愿意写代码的人通常都不太愿意去写文章，由于代码方面的艺术和文字中的美学每每很难兼得，二者都兼得的人一般都已经被西方极乐世界所收罗，我也是只喜欢写代码，让那些字母组成美妙的歌曲，而后自我沉浸在其中自得其乐。而今天，在清明之际，在踏青时节，我仍是忍不住停下来歇歇脚，稍微共享一下最近一直研究的一个很是基础的算法和应用 - 多目标多角度的模板匹配。

       模板匹配，这是一个几十年来一直为业界所重点研究和处理的算法，存在于各类不一样的机器视觉库中，若是哪个没有提供这个功能，那么他将没法获取你们的承认，也就失去了最基本的活力。能够说模板匹配基于机器视觉就至关于数组在编程语言中同样，基础可是不可或缺。

       在2004年时，个人毕业设计中一个很重要的部分也是模板匹配，当时用模板匹配找到每一个量杯中黄色的油的位置，如今看来用那个算法也是醉了，不过能顺利毕业还考得就是他。

       在个人早期博客中，有一篇文章已经谈到了这个算法，详见：标准的基于欧式距离的模板匹配算法优源码化和实现(附源代码)， 可是这个是个很是慢的过程，并且是单目标无旋转的实现，在实际应用中，这个基本没有啥实际的价值。

       在工业应用场合，有着很是普遍使用场景的是多目标多角度的模板匹配（基本无缩放或轻微缩放），这方面实现的比较好的有halcon、海康、康耐视等，国内也有一些小单位有作研究，并且效果不错。在网络上其实也有比较多的文章谈到了多目标模板匹配，基本上都是基于Opencv实现，良心的说也谈到了一些核心技术，可是仍是皮毛，基本都是一带而过，并且实现的效率也基本是没有什么实用价值的，多是怕说多了别人学会了吧。

       虽然在个人实现中，也参考了很多网络上的文章，可是大部分的细节仍是靠的本身的思考和朋友的一些指导，为了尊重他人，我也不打算特别深刻讲解个人实现，可是仍是把一些具备必定深度的问题提出来，也算是回报网络吧。

      一、概述

       这里先提工业界最为经常使用，也是最为基本的模板匹配方式，基于NCC的灰度模板匹配。

      NCC，全称为Normalized Cross correlation，即归一化互相关系数，  在模板匹配中使用的很是很是普遍，也是众多模板匹配方法中很是耀眼的存在，  这个匹配的理论核心基础公式以下：

                                                                         （1）

        该方法也存在于Opencv的matchTemplate中，较之其余的CV提供的匹配方法，该算法对于光照、噪音等等的影响，稳定性更佳，也是halcon等商用软件内嵌的基于像素的模板匹配标准方法。

        他的理论匹配度范围是[-1，1]，为-1时表示2副图像的极性彻底相反（原图和反色后的图），为1则表示两幅图彻底同样。通常咱们在计算NCC的时候都是取的绝对值，所以，一般NCC的取值为[0,1]，值越大，表示两幅图像越类似。

        实际编程实现时，千万不要直接用这个公式，若是你使用，那你离砸电脑已经不远了，请必定要相信我。

        实际中，咱们都用下面的式子来实现编码（不要问我里面的符号的意思，两个图来自不一样的资料，里面的字母也不同，可是要研究的这个的人都应该能看懂）：

                                         （2）

   　这个式子看上去更为复杂，可是实质上和公式（2）和公式（1）就是同一个东西。公式（2）咱们能够把他拆解为7个部分。咱们一1道来。

       ①、这个留到最后在说。

       ②、T表明的是模板，那么②对于固定的模板来讲就是一个定值，在匹配前能够直接计算好，无需担心耗时问题。

       ③、I 表示的搜索图像中的和模板同样大的一个子块，很明显这个累加有多重方法能够快速的实现，好比比较原始的积分图技术，或者个人BoxBlur里的那种更为快速的实现，这一项也是和参数无关的。

       ④、第四项处理方式同 ②项，无需多言。

       ⑤、第五项彻底同第二项，同时四和五项做为一个总体也能够提早计算好，不参与匹配过程的计算。

       ⑥、第六项处理方式同第三项，也无需多言。

       ⑦、第七项彻底同第三项，直接使用。

       前面的分析代表，第二至第七项要买能够做为常量提早计算好，要么就能够经过某种技术实现O(1)的快速计算，那么如今咱们再回过头来在看第①项，他是模板图像和搜索图像同面积区域像素的一个卷积，这个是没法用某种优化技巧去实现和模板大小无关的快速实现的，注定了他就是NCC计算式中最为耗时的部分。

       有人说卷积能够有FFT实现优化，没错，很是赞成您的观点，可是，朋友，FFT虽然其第一个F表明了Fast，可是呢他在傅里叶的世界是快的，在咱们模板匹配的空间内他受到了一种无形的压迫，在工业界仍是没法接受的。

       所以，咱们注意到在本例中，这个卷积其实都是字节类型的计算，对于一个N*M大小的模板图，这个卷积须要N*M次乘法以及N*M-1次加法，因为是整数计算，自己运算速度还算比较快的，并且若是在PC上我曾经说起过有一种使用SIMD指令的提速方法，大概能有10倍的运力提升，核心是使用_mm_madd_epi16（对应PMADDWD这个指令）。

　　其原理以下：

_m128i _mm_madd_epi16(__m128i a, __m128i b)

     Multiplies the 8 signed 16-bit integers from a by the 8 signed 16-bit integers from b Adds the signed 32-bit integer results pairwise and packs the 4 signed 32-bit integer results. 

     该指令能够一次性执行8个16位数据的乘法及4次加法，而咱们只要提早把8为字节数据转换为16为数据就能够了， 一般这能够由_mm_unpacklo_epi8或者_mm_ cvtepu8_epi16实现。

      二、进一步提高

　  有了上一步的分析，也许你就准备开始动工了，千里之行始于足下嘛，可是做为过来人，我奉劝你除非是为了看效果，不然，你仍是等一等吧，下面的更精彩。

      若是咱们按照公式（2）就开始霸王硬上弓，大部分状况下你须要耐心的盯着你的屏幕鼠标在那里转圈转圈转个10来秒吧，而后就能够见证奇迹了，也许做为初学者，你会心动，而那些须要靠这个吃饭的，可能就是心痛了.........

      为了获得更快的搜索速度，一个最容易想到的策略就是图像金字塔，图像金字塔每增长一层，图像的点数和模板的点数救减小4倍（理论数据，实际因为非2的宽度和高度，不必定正好是4倍），那么不考虑Cache等其余的因素，理论上每增长一层金字塔救能够提速16倍，所以，若是咱们构建了一个4层的金字塔，那么在第4层金字塔上的一次完整匹配，其计算的次数和原始的数据相比，就能减小4096倍。

     咱们先以无旋转单目标为例进行简单的说明，当咱们在金字塔最高层进行一次完整的匹配后，咱们能够找一个全局的极值点，这就是在顶层匹配时的最佳匹配位置，此时，咱们能够将顶层匹配的结果映射到金字塔的下一层中，简单的说就是将找到的匹配点坐标的X/Y值乘以2，那么为了保证在下一层中的精度，此时的搜索区域须要适当的增大，好比选择匹配点周围5*5或者7*7的一个区域，而后在这感兴趣的区域中再进行一个局部的匹配计算，此时只须要计算25或者49次小匹配的计算，当计算完毕后，再次提取出这个小区域的极值，做为本层的最终种子点，重复这个过程，直到金字塔的最底层（即原始搜索图像）为止。

　　稍微分析下，假定上述咱们的搜索的局部范围为5*5大小，金字塔取为4层，那么总体的计算量是原始直接计算的多少呢，这样评价：

                式中M和N分别表示图像的宽度和高度。

       通常来讲M和N都是至少以百为单位的，所以上述计算式的结果至关小，速度能够获得极大的提高。

       这里针仅对一个问题进行展开，即金字塔构建时采用何种下采样算法，讨论以下：

       ①、最近邻，这个结果太粗糙，不利于算法稳定性，能够直接Pass掉。

       ②、双线性插值，这个兼顾速度和效果，是个能够考虑的选项。

       ③、三次立方插值，这个东西在图像放大时是个不错的选项，而金字塔得创建是缩小过程。

       ④、兰索斯插值，这个相似于三次立方插值。

       ⑤、采用严格的高斯金字塔采样矩阵，一个5*5的矩阵，以下所示： 

           

 

        对于缩小，其实③④⑤都不是很好，③内部是一个4*4的取样，④是一个8*8的取样，我举一个简单的例子。

                  

   　　　　原图（注意红色框部分的效果）                                              最近邻　　　　　　　　　　 双线性　　　　　　　   　 三次立方　　　　　　　　　兰索斯

      你们可用大一点的屏幕去观察，能够看到红色方框内在原图部分为很是光滑的，而四个插值中，最近邻有所模糊，三次立方和兰索斯在风车叶片边缘出现了锯齿，只有双线性完美的保存了叶片的边缘光滑性。

      那么通过个人实测，一种更好的方式是直接使用2*2均值下采样，也就是使用2*2区域内的全部像素的平均值，2*2均值滤波器有一个很是好的特性，他没有频率响应问题，而较大的滤波器均存在该问题。同时，使用该滤波器还有一个优异的特性，即他能够以很是高效的方式实现。

      当图像的宽度和高度都为2的整数倍时，若是选用双线性插值创建下一层金字塔，此时的双线性就退化为了2*2均值滤波器。

      这里留一个问题你们共同探讨吧：

      问题1： 使用2*2均值滤波器时，对于非偶数的宽度和高度，好比W0 = 101,下一层金字塔的宽度W1究竟是取50，仍是取51呢，若是取51，那么第51个像素如何获取结果（2*2取样会越界）呢？

      三、多目标

      现实中一般须要在一副原图中寻找多个匹配的目标，此时，咱们的难度就有所上升了，金字塔咱们依旧还可使用，可是如何区分多目标呢。

      在单目标时，咱们对最高层金字塔进行了全图匹配后，只须要取最大的匹配值做为候选点，这也是理所固然的，当问题来到多目标时，一个天然的想法就是对匹配值进行排序，而后取前N个最大值做为候选点。可是，这种直接的想法确实很是的不科学的，由于一般在某个最大值附近，因为邻域的类似性，还存在不少和最大值很是接近的得分点，而他们实际上都是对应同一个目标。

      在百度搜索多目标模板匹配，大部分状况你会找到以下的一段代码：

Point getNextMinLoc(Mat &result, Point minLoc, int maxValue, int templatW, int templatH) {
    int startX = minLoc.x - templatW / 3; int startY = minLoc.y - templatH / 3; int endX = minLoc.x + templatW / 3; int endY = minLoc.y + templatH / 3; if (startX < 0 || startY < 0) { startX = 0; startY = 0; } if (endX > result.cols - 1 || endY > result.rows - 1) { endX = result.cols - 1; endY = result.rows - 1; } int y, x; for (y = startY; y < endY; y++) { for (x = startX; x < endX; x++) { float *data = result.ptr<float>(y); data[x] = maxValue; } } double new_minValue, new_maxValue; Point new_minLoc, new_maxLoc; minMaxLoc(result, &new_minValue, &new_maxValue, &new_minLoc, &new_maxLoc);
    return new_minLoc; }

　　这是一段比较好的参考代码，可是仅仅是比较好，还不能解决不少问题，可是咱们能够从中学习到一些东西。

        代码中输入参数中有一个参数是前一次的最小值的坐标，而后在这个坐标附近必定的矩形范围内（上述代码是模板图像的1/3尺寸），将得分值修改成某个很大的值，接着再进行全范围的最大和最小值定位，此时，确定就定位到了离输入最小值有必定距离的另一个最小值，而输入最小值附近的那些次小值就不会干扰结果。

       注意上面代码是最小值，由于他用的检测指标是CV_TM_SQDIFF_NORMED，而非NCC，对于NCC，则须要归为最大值。

       这里的templatW / 3和templatW / 3有点重叠范围的意思了，可是还彻底不够。

       在最顶层金字塔中找到了多个目标的粗糙位置后，就能够和前所述的同样的方式一步一步的向下一层金字塔进行细化，直处处理到顶层金字塔为止。

       四、多目标+多角度

       当问题来到这里时，整个的难度就是阶跃式的提升了一个档次。

       若是目标存在旋转，为了能找到发生旋转的物体，咱们能够建立多个方向的旋转对象，也就是说，将搜索空间离散化，此时，有两个可选的方式：一个是旋转搜索图像，而后用模板在旋转后的图像中搜索，二是旋转模板，用旋转后的模板在搜索图像中定位。咱们说，第一种方式基本不可取，缘由有三。

       (1)、搜索图像通常来讲都是较大的图，对其进行旋转耗时比较可观。、

       (2)、实际状况须要多个角度的旋转，对原图旋转内存方面也会有过多的消耗

       (3)、工业应用时，通常模板比较固定，而搜索图像老是时刻变化的。

       当选择第二种方法时，对于较小的模板图像，是能够在执行搜索前把相关旋转信息提早准备好，在搜索时刻直接使用，而无需作无谓的耗时。

       此时，在金字塔的最顶层，须要作的计算工做也有所增长，咱们须要对每一个角度的模板都作一个全图的匹配，获得匹配的结果，而后对每一个可能点，选择匹配度最大的那个角度做为顶层的候选点。

       相似的，在向下一层金字塔映射时，不只仅须要映射匹配点的X和Y坐标，还须要映射角度信息，同理，为了保证角度方面的精度，也须要适当的扩大角度的搜索区间。

       五、更多的问题

       实际上，为了实现多目标和多角度的匹配，还存在不少不少的细节问题，须要取研究，这些方面的细节有些已经有了部分结论，可是大部分在网络上鲜有探讨和方向，这里列出一些问题和你们共同探讨一番。

       问题2：金字塔多少层比较合适？

       前面提到了金字塔层数越多，所需的计算量就越少，可是同时带来一个问题，就是计算的精度会愈来愈差，这是由于，随着金字塔层数的增长，由于二次采样图像会不断的获得平滑，在图像分辨率变得过低时一些基本的特征已经彻底丢失，各类原本不想关的信息已经彻底融合在一块儿了。所以，一个合适的金字塔层数就尤其重要。一些成熟的商业软件通常可提供用户自行输入金字塔层数，或者自动肯定。对于大部分用户而言，他们不懂得更多的技术细节，自动设置显得尤其重要。

       咱们以工业界最为出色Halcon软件为例，通过屡次测试，他的金字塔层数的自动设置是很是智能化的，基本上自动能够保证速度最优同时效果稳定，一般，咱们认为金字塔的层数只和模板图像的尺寸有关，可是，一些仔细的测试代表，两幅一样大小的模板图，halcon也有可能但不必定返回不一样的值，估计这其中还用到了图像的一些类似性或者结构方面的一些信息做出的综合判断，不过，一个最基本的规律，仍是能够分享下，那就是若是金字塔某层的短边像素小于或者等于8了，大家这确定是此模板金字塔层数的极限了。

       问题3：角度离散化的间距如何设置最为合理？

       相似于金字塔，角度离散化的间距越大，须要计算的旋转方面的信息就越少，速度则能够更快，可是一样，精度就越差。

       通常来讲，若是模板越大，离散化的间距则须要越小，这是由于较大的模板可以区别更小角度的变化。一般，对于大小100像素的模板，离散的角度步幅可设置为1度。一样，如何自动的设置参数也应该是一个成熟软件的标志，一个可行的建议是离散的步长须要保证旋转时两个相邻模板之间的长或宽的最小差别值不小于1个像素。

       问题4：模板的旋转如何处理？

       一般能够用双线性插值或者三次立方插值，来获取旋转后的数据，不建议用最近邻算法，可是不一样的旋转算法，最后获得的匹配结果会有所不一样，同时这也就说呢，其实带角度的模板匹配，理论上很难获取精确解，由于你毕竟不知道原始的旋转算法是何种，好比我从一个未旋转图像中扣一个矩形出来做为模板， 而后把图像旋转各10度，用halcon对模板进行匹配，获得的结果哪怕选择亚像素也不会是精确的10度。

      问题5：旋转后无效数据的部分怎么办？

      在对模板进行旋转时，除非是几个特殊的角度，好比0、90、270、360度，不会产生额外的信息，其余的角度，都会有一些未知的区域存在，以下图所示：

  

         　　　　　　原   图                                                                            旋转必定的角度 （双线性插值）                                                   边缘处局部放大

       中间为旋转后的模板，红色部分为新出现的未知区域，若是说咱们这个图做为一个总体，去和原图像进行匹配，也就是说让红色部分参与了匹配，这很明显会获得一个得分较低的错误结果。

       一种也许可行的方式是，咱们把这些红色的区域剔除在匹配有效的数据以外，这时，又会面临另一个新的问题，在图像的边缘部分如何处理。在上面边缘处局部的方法图中，咱们能够看到，因为插值的特性，边缘处未能在原始图像中采集到足够的采样点，所以，选择了红色背景色做为融合的基色，此时的结果像素就不彻底是属于原始图像了，怎么办？我的的一个建议是既然这部分像素也被污染了，那就不该该参与后续的匹配。

 　  此时，在编程方面就须要继续克服下一个困难了，前面概述方面讲的一些优化方式又不能直接使用了， 真是他妈的痛苦，因此，借用某一本书里的经典劝退名言吧：聪明的机器视觉用户都依赖标准软件包来提供这些功能而不会试着本身实现这些算法。

      问题6：各层金字塔的角度离散值如何分配？

      一般，在金子塔的最底层（和原图同样大小那一层），可按照前述的自动角度幅值来一步一步的旋转图像，而后随着金字塔的层数增长，根据模板在每层金字塔中都会缩小2倍的这个事实，在相应金字塔上模板的角度步幅也能够增长2倍，所以，若是在金字塔最底层上使用的角度步幅喂1度，那么在第四层上可使用8度做为角度步幅。

      这样作带来的好处有不少，由于，一般，咱们须要在最顶层作多角度的全图的匹配，这个计算的相对来讲比较大，若是角度步幅在该层得以以指数级别减小，则计算量也会有量级的下降。

      问题7：各层金字塔的最低得分值如何肯定？

      前面一直没提这个得分的问题，在单目标中通常是不存在这个问题的，当有多个目标时，由于可能不是完整的匹配，所以通常须要客户肯定一个最小的得分值，固然这个得分是指的在最底层时的值。当有了金字塔后，由于下采样的一些因素，为了有更高的容错率，通常是建议每增长一层金字塔，最小的得分值须要适当的下降。好比，用户设置的最小得分为0.8，后续各层的得分能够为：

       　　　　 0.8          　　0.8*0.9　　　　0.8*0.9*0.9　　　　　0.8*0.9*0.9 * 0.9　　

　　固然，其余的调整方式也何尝不能够，可是，整体的区域需保证最小得分愈来愈低。　　

　　问题8：MaxOverlap是什么鬼，内部是如何操做的？

       在Halcon中，还有很是重要的参数-MaxOverlap，一个介于0和1之间的参数，前面也一直没有谈到过。其实，在真正的应用中，存在着一些目标之间有必定程度重叠的状况，这个时候，若是按照前面的那些处理方式，通常就只能获取到重叠对象的某一个，而丢掉了其余的信息。当有了MaxOverlap参数时，咱们就能够根据这些对象的重叠重复，来决定是否剔除掉某一个不须要的目标。

       可是，这也是一个比较难以琢磨的对象，Halcon帮助文档中的说法是取的对消的最小外接矩形中的重叠比例。说实在的，编程时，这个规则应该还不够。好比，若是有4个对象，他们都有所重叠，请问这个时候，这个MaxOverlap是指的全部的重叠量合并在一块儿呢，仍是最大的重叠，抑或是按照得分顺序第一个和之重叠的呢。目前，我也是对这个比较摇摆。

      同时，另一个比较难以肯定问题是，这个重叠是在金字塔的最顶层就进行判断仍是如何呢，若是在金字塔的最顶层不进行判断，那么金字塔顶层中得分大于MinScore的则有不少个点，若是把这些点都直接向上一层进行映射，那个计算量也是至关客观的。

      问题9：亚像素坐标和角度是一块儿执行的吗，仍是分开的？

      没有亚像素的模板匹配是没有灵魂的，特别是带有角度的匹配。由于，正如前面所述，咱们对角度采用了离散化。那么这个时候计算的最终角度结果必然是离散化后的序列里的某个值，这样的精度有的时候是不够的。

      经过亚像素技术，寻找局部区域里的最大值，课得到更高的位置精度和角度精度，则能够有效的得到精准的定位。

      可是目前仍是不清楚坐标+角度组成的4维的空间的亚像素计算公式是如何的。

     问题10：速度优化？

      以上谈的一些优化基本上是结构上的优化，或者说原理上的优化，在编码上还能够考虑时候用SIMD指令进行处理，尤为是全图匹配的计算过程。

      六、公关成果

       大概前先后后再一块儿有折腾了大概一年的时间吧，初步仍是有了必定的成果，不管是在速度上，仍是在准确度中，仍是能解决必定的工程问题了。

       共享一些测试图，和你们一块儿比较。

                   

     　　模板图　　　　　　　　　　　　　　　　　待搜索图　　　　　　　　　　　　　　　　　　　　　　　　　　        搜索结果

 

                  　　

         模板图　　　　　　　　　　　　　　　 　　待搜索图　　　　　　　　　　　　　　　　　　　　　　　　　　        搜索结果

      基于NCC的速度方面和不少参数有关，好比角度的搜索范围，金字塔层数，模板的大小（通常模板大，速度反而快些，特别小的模板则很是耗时，知道缘由吗？），重叠的大小等等都有关。一些简单的测试数据以下：

        　　

       虽然历经千辛万苦，在磨砺了好久以后，也对这个初有小成，基本实现了这样的一些过程，可是和halcon相比，不管是从稳定性仍是效率方面都仍是有必定的差距的，因此标题中的无限接近 就是一句诳语而已。

      本算法目前已经集成到国产视觉软件Malcon中，详情请看 中国的Malcon跟德国的Halcon的相比的优缺点 。

      Malcon官方博客：https://www.cnblogs.com/Malcon

      或者点击：     　   https://blog.csdn.net/lindrs/article/details/114113280?spm=1001.2014.3001.5502

      Malcon中还集成了不少其余经常使用的200多个算子，相信您一种能从其中找到你所须要的。

      若是你但愿有一个简单的可视化测试界面，能够从以下连接中获取，可是请注意这个Demo自己是有一些BUG的（不影响测试使用），请不要将其直接应用到工业环境中，以避免形成没必要要的损失。

       可视化测试Demo:   https://files.cnblogs.com/files/Imageshop/TemplateMatching.rar

       七、更多疑惑

       弄得越多，发现不了解的也越多，特别是在研究比对Halcon的过程当中，还发现了一些暂时没有弄清楚的事情，好比：

       ①、在不勾选亚像素时，Halcon返回的坐标值不少状况下也是带小数点（特别是非0度时的结果），这个做何解释。

       ②、当模板的理论角度是0度时，即便按照AngleStart和AngleStep的值依次计算，取样的角度值不会准确的取为0，并且也未勾选亚像素，他也能正确的返回0值，难道他对0度作了特别的处理吗？

       这些还须要慢慢的取探索。

       若是您以为本博文对您有所帮助，也可慷慨解囊。