Linemod 代码笔记

最近了解到 Linemod 这个模板匹配算法，印象不错
准备仔细学习一下，先做点代码笔记，免得后面不好回顾
目前的笔记基本上把 核心流程都分析得比较清楚了，除了一些阈值的选取

opencv 的contrib 模块有这个算法的实现

我看的代码来自这里
https://github.com/meiqua/shape_based_matching

先大概记录下 代码思路：
分两个阶段， train 和 test

Train

Train 中 ， shapeInfo_producer 负责用来对 模板进行 各种旋转和尺度缩放，
shapes.src_of 可以根据旋转和尺度 生成变换后的 模板

对每一个模板 执行 detector.addTemplate 操作，

最后调用 shapes.save_infos 和 detector.writeClasses 这两个保存训练 结果。保存的信息用于 后续的匹配中。

首先构造
line2Dup::Detector detector(20, { 4, 8 });
第一个参数为 特征点个数 ， 第二个参数是一个 vector, 每个元素代表每一层的T
构建 this->modality 对象

shape_based_matching::shapeInfo_producer shapes(padded_img, padded_mask);
两个入参都是 图像，第一个是用 输入图像构建，填充像素为0， 第二个用输入图像大小的大小构建掩码图像，掩码为1， 填充像素为0

然后填充shapes.scale_range、 shapes.scale_step、 shapes.angle_range 、shapes.angle_step
这四个是对模板图像进行 尺度缩放 和 旋转的 量

shapes.produce_infos();
主要是用 尺度范围 和 旋转范围 的组合 构建 std::vector infos
然后 就是 遍历 shapes.infos
执行
detector.addTemplate(shapes.src_of(info), class_id, shapes.mask_of(info));

shapes.src_of(info) 产生变换后的图像
class_id 是一个固定的字符串
shapes.mask_of(info) 返回 shapes.src_of(info) 产生变换后的图像是否大于0的 掩码图像
addTemplate 是 核心函数，主要作用为 提取模板图像的特征点，即梯度较强的点，得到 这些点的坐标和梯度方向值。

接着调用两个函数

1. shapes.save_infos 保存 的信息是 每张图片是 原始图像经过哪种旋转和缩放得到的
2. detector.writeClasses 则 保存 每个模板 的信息，包括cropTemplates(tp) 后的高宽和坐标、 特征点坐标信息，特征点的label 就是梯度方向

=============================================================================

Detector::addTemplate

1 modality->process(source, object_mask)

这个是 直接构造一个 ColorGradientPyramid 对象，返回其指针
ColorGradientPyramid 构造函数中 update(); ，内部是
quantizedOrientations(src, magnitude, angle, weak_threshold);
先做 高斯模糊， 然后 在水平和垂直方向 调用 Sobel，
调用 phase 计算梯度方向，

调用 hysteresisGradient， 主要输出就是 quantized_angle
过程为： 先把 连续的梯度方向 划分为16个区间， 然后量化为8个方向
quant_r[c] &= 7; 这个代码还没看明白，这 相当于把一个整数 对8 求模
这么做没问题应该是因为 认为 180度和190度之间的方向 和0度到10度之间的 方向是一个方向。

然后就是 对梯度幅值 超过一定阈值的 像素点 的 3*3 邻域 求 梯度直方图
投票数 超过 阈值的 方向 作为最终的 量化方向

至此， modality->process 完成
返回 一个 Ptr qp

然后 开始遍历金字塔每一层， 如果不是最底层， 那么 qp 降采样，并且 做梯度量化操作， 即调用上面的 update()

然后qp->extractTemplate(tp[l])

这一步是 提取第 L 层特征点， 保存在 tp[l]中。 细节参考后文
说明： tp是个vector， 每个 元素都是一个模板，对应金字塔某一层提取出来的特征点

每一层都遍历完后， cropTemplates(tp)

这个函数 先 遍历每一个 模板， 找出特征点最大最小坐标，注意，高层次的金字塔图像的坐标会进行放大（根据层次）
得到 4个最小、最大坐标。 注意： 是所有层共用信息

然后再一次遍历每个模板， 调整 templ.width ，templ.height ，templ.tl_x，templ.tl_y
然后用 templ.tl_x，templ.tl_y 修正了特征点坐标，
TODO: 这就 有点麻烦了， 修正后的 坐标肯定和 原始图像 对应不上了啊！

返回 Rect(min_x, min_y, max_x - min_x, max_y - min_y)
但 外部并未接收 这个返回值

addTemplate 的最后 template_pyramids.push_back(tp);
ColorGradientPyramid::extractTemplate(Template &templ)
函数输出应该是 templ.features， 即提取出 特征点
先对 mask 进行 腐蚀，

Magnitude 是 之前 quantizedOrientations 中计算出的梯度幅值（梯度平方和）

对 Magnitude 搞一个 遍历，
如果对每个像素，如果 magnitude_valid 值 大于0
如果其邻域内 有像素的梯度幅值超过它，
那么 is_max 为 false, 如果遍历完后 ， is_max 为true, 那么 所有 邻域像素对应 magnitude_valid 值 置为0

通过上述检验的点 ， 如果 幅值超过阈值， 且 方向不为 0， 进入 candidates
（注意 opencv在这里的实现方法， 先设置了一个 score = 0， 如果没通过上述检验， 该值依然为0， 这种实现方法好吗？）

遍历完后，如果 candidates 个数低于阈值， 返回 false, 此次 抽取失败。。。

对 candidates 按照 score 进行一次稳定排序
selectScatteredFeatures 最后 从 candidates 中 选取一些 散得 比较开的点， 这里while 循环写得还比较有技巧， 如果遍历完一轮， 数量不够，那么 降低 距离阈值， 再选！
和 orb-slam或者说opencv 里面 ORBextractor 提取特征点 那个 四叉树的方法谁优谁劣？

选取的特征点保存 在 templ.features 中

Test

先读取 train 阶段保存的两个信息文件
detector.readClasses(ids, prefix + “myCase/%s_templ.yaml”);
读取 每个模板 的信息，包括cropTemplates(tp) 后的高宽和坐标、 特征点坐标信息，特征点的label 就是梯度方向。
构建出： class_templates

shape_based_matching::shapeInfo_producer::load_infos
每张图片是 原始图像经过哪种旋转和缩放得到的

对测试图像 进行一下调整， 使得高宽都是 16 的倍数

auto matches = detector.match(img, 90, ids);
90 是阈值， ids 是 训练时 指定的id字符串 test

然后 modality->process(source, mask)，
这个调用在前面已经介绍过了，会 构造一个 ColorGradientPyramid 对象，对source图像计算量化后的梯度信息

然后遍历 金字塔， construct response map
先不看 具体的函数调用实现过层， 从函数名字 和 注释来看， 这就是 论文当中第三节讲的东西， 包括 方向扩散spread、 梯度响应计算computeResponseMaps、 线性化存储linearize。 最终存在在 LinearMemoryPyramid 结构里面。

遍历class_ids， 从 class_templates获取 对应 std::vector
matchClass(lm_pyramid, sizes, threshold, matches, it->first, it->second);
这个函数完成整个匹配过程

=============================================================================

Detector::matchClass

遍历template_pyramids， 提取出 每个 Template，
调用 similarity， 计算相似性， similarity中， 核心调用是 accessLinearMemory，
这里面第一行代码
const Mat &memory_grid = linear_memories[f.label];
很关键，这是根据模板中特征点 来 定位 response map 相应的数据
定位到以后，然后 就是 SIMD 指令 来 累加数据了！

static void spread(const Mat &src, Mat &dst, int T)

这个地方实现的是 论文3.3 节的所谓 梯度方向展开
所要实现的功能很好理解， 即把每个像素及其邻域的离散化的梯度方向进行 或运算。
OpenCV 这里再一次展现了实现技巧， 最直观的方法是 每次遍历一个像素时，取出其所有邻域内的像素的梯度方向值，然后做一个或运算， 这样做 内存访问性能较低， 因为图像的下一行和上一行 距离较大， 很可能缓存命中失败。

OpenCV 的做法是： 每次遍历时， 只做整个邻域内某个特定位置的像素梯度方向值 的 或运算，这个地方说的邻域包含像素自身，即邻域中心。 所以总共循环 T*T次。 T 为邻域直径。
这样做， 内存访问友好，并且方便使用 SSE指令进行优化， 因为连续参与运算的数据在内存中是连续的！

static void computeResponseMaps

(const Mat &src, std::vector &response_maps)

实现论文3.4节 响应图的计算
这个地方 把论文中的相似度 也给离散化了。
并且事先计算了 某个方向 和 某组方向的余弦值的最大值，并且离散化， （或者称为根据余弦值 实行打分制） 存储到一个数组SIMILARITY_LUT 中，即查找表。 这个查找表中针对某个方向的值有32个元素， 总共8个方向， 所以有 256个元素。 32个元素中 ， 又分为两组， 前16个是8个方向中前4个方向的各种组合 与 当前32个元素针对的方向 的余弦值的最大值对应的得分。

这个数组， 上交这个学生 对原来的值 进行了修改： 1,2–>0 3–>1
为什么这么改？
https://zhuanlan.zhihu.com/p/35683990
这篇文章给出了 修改的解释

论文3.4 节 也给出了 这个查找表的计算啊！

疑问待定： n0 为8的时候， 针对某个方向的查找表元素 按照论文实际上应该是有 2的8次方， 即 256种情况。 这个地方是不想搞出那么大一个数组， 所以， 把8位分拆成两组， 每组只需16个元素， 然后再进行一次比较，拿到最终的最大值？ 为啥不直接构建大小为 256*8的查找表? 这样可以省掉一次 max的运算。
看了下 _mm_shuffle_epi8 的介绍

这个地方 index 只用低4位进行运算， 也就是只支持 4个bit作为索引值，
如果只能用这个指令，的确 只能把 8位拆分成两组4位，再max
不知道有没有 能直接用8位作为 所以索引的SSE指令

static void linearize

(const Mat &response_map, Mat &linearized, int T)

这个是改变存储方式，先行后列， 间隔T 读取，然后写入。没有比较复杂和特殊的处理。

similarity_64

这个函数计算 模板和 输入图像的 相似性， 即论文中的 similarity map
计算相似性的时候， 并不是 把 模板上的每个像素都和 输入图像上对应的像素 一一对应，然后进行 某种计算， 这和 NCC, SSD 这些方法的做法不一样！一开始受这些方法先入为主的影响，导致论文里的Fig 7 以及代码中的操作

实际上， 只比较模板上提取的特征点， 以及 模板 覆盖在 输入图像上某个位置时， 这些模板特征点对应到 输入图像上的像素点 之间的梯度差异。

意识到这点以后，就比较好理解代码了。 因为模板需要在输入图像上进行 滑动，所以产生了 similarity map。 每次滑动，模板和输入图像产生一个 相似度。 模板在 水平和垂直方向进行滑动， 所以 产生一个 二维的相似度矩阵。这个矩阵的宽 自然就是 输入图像的宽减去模板的宽， 也就是代码中的span_x。 高的情况类似。

代码当中用 template_positions 表示 模板的当前滑动位置。

计算similarity map最直观的方法是：对每个模板位置， 找出所有特征点在输入图像上对应的像素， 计算所有梯度方向的相似性，累加。 然后 处理下一个模板位置。

但代码中的做法是： 对每个特征点，计算出所有模板位置上 这个特征点 和 所有输入图像上对应点的 梯度方向相似性，保存到similarity map中。 然后 计算下一个特征点的相似性，累加到 similarity map中。

整个算法中 不是第一次使用这种思路了。