[图像算法] Linemod算法研究
Linemod算法研究
先了解一下大致工作流程：
ref：Gradient Response Maps for Real-TimeDetection of Textureless Objects
注册过程
注册过程需要提取特征点，后续滑窗的时候就值拿这些特征点在场景图上进行滑窗，而不是传统意义上的图片滑窗，这样可以大大加速匹配过程。
提取特征点可以根据模板图本身的梯度情况，将一些梯度较大的点保留作为模板点，因为图像梯度往往包含着形状相关的信息，并且模板点在基于形状的模板匹配的实际场景中往往二值图居多，基本反应的是形状边缘。
现在考虑匹配阶段。
分数计算
由于后面需要进行滑窗算分，算分的过程可以理解为相似度匹配的过程，那么特征点的信息一方面需要包含位置信息，另一方面也需要包含某种值信息，用于和场景图进行相似度运算。
linemod选择保留的值信息是梯度的方向，因此后续进行相似运算也是判断梯度方向的一致性。
再次的，为了加速运算，梯度方向被进行了量化处理，这样做的好处是将方向离散化，离散有限的数据就可以制表，因此相似运算可以变成查表运算。
如果不考虑离散化，梯度的相似度可以理解为两个向量的相似度，即\(cos(\theta_1 - \theta_2)\)，但是考虑到可能真实场景中会遇到一些黑白invert的情形，因此对这一结果加了绝对值，当然了，这种方式并非绝对的好。
离散化之后，可以人为对一些情况进行赋分：
在一些code中，红色方向和红色方向的相似度赋值为4分，红色和蓝色为3分，以此类推，对于大于九十度的，根据公式也是有分的，毕竟是cos绝对值嘛，opencv官方code是给了分的，而在某些code实现中，大于90度的是认为0分。
具体的：
OpenCV：
meiqua实现：
我个人是支持90度以上0分处理的，毕竟invert并不是对称，0分合理。
关于制表的过程，后面细讲。
此外，滑窗匹配的过程中你并不能太绝对的对比特征点和其对应的点的相似度，宽松化处理是要比对特征点和其对应点周围的点的相似度，这样能够增加命中率。
梯度扩散
不考虑梯度扩散，我们直接用naive的滑窗匹配其实就已经结束了，梯度扩散的目的实际上是为了减少漏匹配，在将某一位置的梯度结果扩散给相邻位置，这样在后续匹配的时候其他位置也会有一定可能被匹配到，因此可以理解为宽松化的处理。
需要注意的是，仅对场景图做扩散处理，因为我们在匹配阶段模板就只有特征点这一信息了，到这一节点已经没有模板的图像信息了。
这一过程可以形象的用论文中的图表示：
即原本没有任何梯度信息的一些点也被赋予了梯度信息。
离散化处理
如上所说，离散化处理是为了查表加速。
比较简单，以5位二进制表示为例：
因此上述扩散图的二进制表示就变成了：
问题1：
咱们观察一下该二进制表示的特点，有些角度是组合的比如00011，按照非离散的理解，这俩方向应该可以合成一个方向表示，如果我们合成一个方向表示，会有问题吗？换言之，我们在梯度量化前进行梯度扩散，然后进行量化，理论上就没有这些组合情况了，在后续查表的时候就不需要考虑组合情况的制表，理论运算应该少很多呀。
问题2：
实际spread的过程是从哪里开始spread的呢，随机，还是直接选择左上角点？事实上spread起始点的选取对结果是有影响的，这里需要一定设计。
先解答一下问题1，二进制表示后进行梯度扩散的底层实现是位或实现的，因此速度快。
对于问题2，code中都选择的是从右下角开始扩散，扩散到左上角，这样信息会向坐上偏移，可能更符合对坐标的量化，即向下取整，实际实现的时候代码也是很高效的，后续探讨。
提前制表
前面说了，如果场景图完成了量化，那么理论上你的特征点在任意位置的得分都是可以穷举的，按照8位量化讨论，你也就8个方向，我提前算好场景图每个位置的8个方向的相似度，到时候你模板滑到这了，我就拿你模板的方向去我的表里查，立马就知道结果了。
因此，我们只需要8个2维表就可以完成这一任务了。
事实上，如果先完成场景图的量化再制表，就失去意义了，因为相当于你在匹配的时候花了更长的时间来做表格，之后再用。。。
所以，制表是发生在你写代码阶段的。
场景图的情况显然就不止8种方向了，因为场景图会存在因为spread过程导致的梯度叠加，8位的话就有\(2^8\)种叠加情况，因此你需要制作一个\(2^8 * 8 = 2048\)大小的表格。
但是实际上我们在看code的时候发现，人家8位用的是256的表格，这是一个优化。
在制作表格的时候，分为高8位和低八位分别制表，在匹配的时候可以低8位查一个分，高八位查一个分，最终结果取两者的最大值，和直接查2048的表格结果是一样的，更加形象的理解可以理解为向量可以分解为左半区间和右半区间两个分量，由于只考虑方向，所以你取两个分量分别相似度的最值即可。（当然了，这种理解其实经不起推敲。。）
因此最终制作的表格大小就是\(2^4 * 8 * 2 =256\)
线性化内存
先考虑一下什么时候需要线性化内存。
现代cpu读取数据的方式是一次加载他之后的若干数据位到cache中，称之为cache line，所以如果你要处理的数据不在一个cache line中，就意味着要加载多次数据，会造成浪费，称之为cache miss，因此一般通过内存重排，将较近可能处理到的数据放在一个cache line内，可以大大减少访存次数。
由于linemod要实现对所有点的滑窗，因此实现过程中可以一个点一个点的滑窗，所以数据间隔就是一个滑窗的stride，所以把每隔stride的数据排在一起，这样就能得到 stride x stride个重排数据了，每次丢一个重排数据进去处理即可。
实际的工程化过程中有一些优化如指令集加速等。
Coding
先来一个easist code，仅仅几十行。写了比较详细的注释，适合学习使用，此外，郑重说一下，这个code没有任何加速和特殊处理，不适合工业场景种使用哈。。
这样咱们在对比了加速方法之后，才能非常明显的感受到加速有多明显。。。
#include <opencv2/opencv.hpp>
#include <iostream>
static const unsigned char SIMILARITY_LUT[256] = {
0, 4, 3, 4, 2, 4, 3, 4, 1, 4, 3, 4, 2, 4, 3, 4, 0, 0, 1, 1, 2, 2, 2, 2,
3, 3, 3, 3, 3, 3, 3, 3, 0, 3, 4, 4, 3, 3, 4, 4, 2, 3, 4, 4, 3, 3, 4, 4,
0, 1, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 3, 3, 4, 4, 4, 4,
3, 3, 3, 3, 4, 4, 4, 4, 0, 2, 1, 2, 0, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 0, 3, 2, 3, 1, 3, 2, 3,
0, 3, 2, 3, 1, 3, 2, 3, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
0, 4, 3, 4, 2, 4, 3, 4, 1, 4, 3, 4, 2, 4, 3, 4, 0, 1, 0, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 2, 2, 0, 3, 4, 4, 3, 3, 4, 4, 2, 3, 4, 4, 3, 3, 4, 4,
0, 2, 1, 2, 0, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 0, 2, 3, 3, 4, 4, 4, 4,
3, 3, 3, 3, 4, 4, 4, 4, 0, 3, 2, 3, 1, 3, 2, 3, 0, 3, 2, 3, 1, 3, 2, 3,
0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4};
void Quantize(cv::Mat& src, cv::Mat& dst, std::vector<std::pair<cv::Point, int>>& pts){
cv::Mat smoothed;
cv::GaussianBlur(src, smoothed, cv::Size(7, 7), 0, 0,
cv::BORDER_REPLICATE);
cv::Mat sobel_dx, sobel_dy, angle, angle_q;
cv::Sobel(smoothed, sobel_dx, CV_32F, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
cv::Sobel(smoothed, sobel_dy, CV_32F, 0, 1, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
cv::Mat magnitude = sobel_dx.mul(sobel_dx) + sobel_dy.mul(sobel_dy);
cv::phase(sobel_dx, sobel_dy, angle, true);
// 到范围0-16
angle.convertTo(angle_q, CV_8U, 16 / 360.f);
cv::Size s = angle_q.size();
// 把180度以上的置0, 梯度太小的去掉
for(int Y=0;Y<s.height;++Y){
for(int X=0;X<s.width;++X){
if(magnitude.at<float>(Y,X) < 100.f){
angle_q.at<uchar>(Y,X)=0;
}else{
uchar label = angle_q.at<uchar>(Y,X)&7;
angle_q.at<uchar>(Y,X) = 1 << label; // 到0-8
// label = 0 对应 0000 0001
// label = 1 对应 0000 0010
// ...
pts.push_back({cv::Point(X,Y), int(label)});
}
}
}
dst = angle_q;
}
void Spread(cv::Mat& src, cv::Mat& dst, int T){
cv::Size s = src.size();
dst = src.clone();
for(int TX=0;TX<T;TX++){
for(int TY=0;TY<T;TY++){
for(int Y=0;Y<s.height - TY;Y++){
for(int X=0;X<s.width-TX;X++){
dst.at<uchar>(Y,X) |= src.at<uchar>(Y + TY, X + TX);
}
}
}
}
}
void Detect(cv::Mat& scene_mat,
cv::Mat& temp_mat){
cv::Mat temp_q, scene_q, scene_q_spread;
std::vector<std::pair<cv::Point, int>> feature_pts, tmp;
// 注册模板 -> 提取特征点
Quantize(temp_mat, temp_q, feature_pts);
// 滑窗检测
Quantize(scene_mat, scene_q, tmp);
Spread(scene_q, scene_q_spread, 1);
std::vector<cv::Mat> response_maps(8);
// 生成cvmat格式response_maps
for(int i=0;i<8;++i){
cv::Size s = scene_q.size();
cv::Mat response_map(s, CV_8U);
for(int Y=0;Y<s.height;++Y){
for(int X=0;X<s.width;++X){
uchar val = scene_q_spread.at<uchar>(Y,X);
// 低4位16种组合，高4位16种组合
int idx_low = i * (16+16) + (val & 15); // 0000 1111
int idx_high = i * (16 + 16) + 16 + ((val & 240) >> 4); // 1111 0000
// 取最大
response_map.at<uchar>(Y,X) = std::max(SIMILARITY_LUT[idx_low], SIMILARITY_LUT[idx_high]);
}
}
response_maps[i] = response_map;
}
// 检测
int res_w = scene_mat.cols - temp_mat.cols;
int res_h = scene_mat.rows - temp_mat.rows;
cv::Mat res(cv::Size(res_w, res_h), CV_32F);
for(int Y=0;Y<res_h;++Y){
for(int X=0;X<res_w;++X){
float val = 0;
for(auto& feat:feature_pts){
cv::Point p = feat.first;
int label = feat.second;
cv::Mat& response_map = response_maps[label];
val += float(response_map.at<uchar>(Y + p.y, X + p.x));
}
res.at<float>(Y,X) = val;
}
}
cv::Mat m;
double minVal, maxVal;
cv::Point minLoc, maxLoc;
cv::minMaxLoc(res, &minVal, &maxVal, &minLoc, &maxLoc);
res.convertTo(m, CV_16U, 65535.f / maxVal);
cv::imwrite("C:/Users/xueaoru/Desktop/xx.png", m);
int num_features = feature_pts.size();
std::cout << "Loc: " << maxLoc.x << "," << maxLoc.y <<" Score:" << 100 * maxVal / (4.f * num_features) << std::endl;
cv::rectangle(scene_mat, cv::Rect(maxLoc.x, maxLoc.y, temp_mat.cols, temp_mat.rows),cv::Scalar(128), 1);
cv::imshow("res", scene_mat);
cv::waitKey();
}
int main(){
cv::Mat temp_mat = cv::imread("C:/Users/xueaoru/Desktop/test_job/temp_image7.png", 0);
cv::Mat scene_mat = cv::imread("C:/Users/xueaoru/Desktop/test_job/gray7.png", 0);
cv::resize(temp_mat, temp_mat, cv::Size(), 0.25, 0.25);
cv::resize(scene_mat, scene_mat, cv::Size(), 0.25, 0.25);
Detect(scene_mat, temp_mat);
return 0;
}
大致效果(最后一行为response_map)：
Improvements
主要讲解下meiqua的code中一些加速部分。
- 金字塔加速
大致思路是先用一个下采样的场景图进行匹配，得到若干proposal，再在原始尺度的proposal附近再次匹配，这样的好处是低level的图像处理快，但是不太准确，原始尺度只需要在proposal附近重新寻找就好了，速度上是有明显提升的。
目前在测试中只是发现场景图分辨率太小时下采样的图由于spread等因素会有点小问题。。
void Detector::matchClass(const LinearMemoryPyramid &lm_pyramid,
const std::vector<Size> &sizes,
float threshold, std::vector<Match> &matches,
const std::string &class_id,
const std::vector<TemplatePyramid> &template_pyramids) const
{
#pragma omp declare reduction \
(omp_insert: std::vector<Match>: omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))
#pragma omp parallel for reduction(omp_insert:matches)
for (size_t template_id = 0; template_id < template_pyramids.size(); ++template_id)
{
const TemplatePyramid &tp = template_pyramids[template_id];
// First match over the whole image at the lowest pyramid level
/// @todo Factor this out into separate function
const std::vector<LinearMemories> &lowest_lm = lm_pyramid.back();
std::vector<Match> candidates;
{
// Compute similarity maps for each ColorGradient at lowest pyramid level
Mat similarities;
int lowest_start = static_cast<int>(tp.size() - 1);
int lowest_T = T_at_level.back();
int num_features = 0;
{
const Template &templ = tp[lowest_start];
num_features += static_cast<int>(templ.features.size());
if (templ.features.size() < 64){
similarity_64(lowest_lm[0], templ, similarities, sizes.back(), lowest_T);
similarities.convertTo(similarities, CV_16U);
}else if (templ.features.size() < 8192){
similarity(lowest_lm[0], templ, similarities, sizes.back(), lowest_T);
}else{
CV_Error(Error::StsBadArg, "feature size too large");
}
}
// Find initial matches
for (int r = 0; r < similarities.rows; ++r)
{
ushort *row = similarities.ptr<ushort>(r);
for (int c = 0; c < similarities.cols; ++c)
{
int raw_score = row[c];
float score = (raw_score * 100.f) / (4 * num_features);
if (score > threshold)
{
int offset = lowest_T / 2 + (lowest_T % 2 - 1);
int x = c * lowest_T + offset;
int y = r * lowest_T + offset;
candidates.push_back(Match(x, y, score, class_id, static_cast<int>(template_id)));
}
}
}
}
// Locally refine each match by marching up the pyramid
for (int l = pyramid_levels - 2; l >= 0; --l)
{
const std::vector<LinearMemories> &lms = lm_pyramid[l];
int T = T_at_level[l];
int start = static_cast<int>(l);
Size size = sizes[l];
int border = 8 * T;
int offset = T / 2 + (T % 2 - 1);
int max_x = size.width - tp[start].width - border;
int max_y = size.height - tp[start].height - border;
Mat similarities2;
for (int m = 0; m < (int)candidates.size(); ++m)
{
Match &match2 = candidates[m];
int x = match2.x * 2 + 1; /// @todo Support other pyramid distance
int y = match2.y * 2 + 1;
// Require 8 (reduced) row/cols to the up/left
x = std::max(x, border);
y = std::max(y, border);
// Require 8 (reduced) row/cols to the down/left, plus the template size
x = std::min(x, max_x);
y = std::min(y, max_y);
// Compute local similarity maps for each ColorGradient
int numFeatures = 0;
{
const Template &templ = tp[start];
numFeatures += static_cast<int>(templ.features.size());
if (templ.features.size() < 64){
similarityLocal_64(lms[0], templ, similarities2, size, T, Point(x, y));
similarities2.convertTo(similarities2, CV_16U);
}else if (templ.features.size() < 8192){
similarityLocal(lms[0], templ, similarities2, size, T, Point(x, y));
}else{
CV_Error(Error::StsBadArg, "feature size too large");
}
}
// Find best local adjustment
float best_score = 0;
int best_r = -1, best_c = -1;
for (int r = 0; r < similarities2.rows; ++r)
{
ushort *row = similarities2.ptr<ushort>(r);
for (int c = 0; c < similarities2.cols; ++c)
{
int score_int = row[c];
float score = (score_int * 100.f) / (4 * numFeatures);
if (score > best_score)
{
best_score = score;
best_r = r;
best_c = c;
}
}
}
// Update current match
match2.similarity = best_score;
match2.x = (x / T - 8 + best_c) * T + offset;
match2.y = (y / T - 8 + best_r) * T + offset;
}
// Filter out any matches that drop below the similarity threshold
std::vector<Match>::iterator new_end = std::remove_if(candidates.begin(), candidates.end(),
MatchPredicate(threshold));
candidates.erase(new_end, candidates.end());
}
matches.insert(matches.end(), candidates.begin(), candidates.end());
}
}
- 线性内存
按照每stride T选择一个元素进行重排，最终可以得到T x T个重排后的内存片段，之后每个内存片段会单独使用，具体是用MIPP进行指令集加速，每次同时加载N个数据到寄存器进行处理，具体见相似度计算的代码。
static void linearize(const Mat &response_map, Mat &linearized, int T)
{
CV_Assert(response_map.rows % T == 0);
CV_Assert(response_map.cols % T == 0);
// linearized has T^2 rows, where each row is a linear memory
int mem_width = response_map.cols / T;
int mem_height = response_map.rows / T;
linearized.create(T * T, mem_width * mem_height, CV_8U);
// Outer two for loops iterate over top-left T^2 starting pixels
int index = 0;
for (int r_start = 0; r_start < T; ++r_start)
{
for (int c_start = 0; c_start < T; ++c_start)
{
uchar *memory = linearized.ptr(index);
++index;
// Inner two loops copy every T-th pixel into the linear memory
for (int r = r_start; r < response_map.rows; r += T)
{
const uchar *response_data = response_map.ptr(r);
for (int c = c_start; c < response_map.cols; c += T)
*memory++ = response_data[c];
}
}
}
}
这里相似度运算其实是有重复运算，结果是每个特征点和所有场景图上点的相似度，但右边有些地方模版图是超出场景图的，所以后续会用border给限制住。之所以这里重复运算，也是为了指令集加速，这就是为啥非得每个特征进行滑窗，而不是我的code中对整个模版滑窗那样。。
这里的mipp其实有个疑惑的地方，在前面spread or操作的时候做了内存对齐，而这里没做内存对齐。而实际上前面spread的时候是针对每行行首做的，而前面对图像的pad其实做了行对齐，而opencv在申请内存的时候就有对齐，这就导致了前面spread的内存对齐实际没作用，我debug的时候也发现确实没往那个while循环里走。而similarity函数这里其实应该在mipp前面对齐一下的。。。
此外，内存对齐这里也并不是固定16对齐的，取决于你采用的cpu指令集，如果是128位向量寄存器，那么显然128 / 8 = 16是16对齐，而在我的机器上我的向量寄存器有256位，因此应该是32对齐。
static void similarity(const std::vector<Mat> &linear_memories, const Template &templ,
Mat &dst, Size size, int T)
{
// we only have one modality, so 8192*2, due to mipp, back to 8192
CV_Assert(templ.features.size() < 8192);
// Decimate input image size by factor of T
int W = size.width / T;
int H = size.height / T;
// Feature dimensions, decimated by factor T and rounded up
int wf = (templ.width - 1) / T + 1;
int hf = (templ.height - 1) / T + 1;
// Span is the range over which we can shift the template around the input image
int span_x = W - wf;
int span_y = H - hf;
int template_positions = span_y * W + span_x + 1; // why add 1?
dst = Mat::zeros(H, W, CV_16U);
short *dst_ptr = dst.ptr<short>();
mipp::Reg<uint8_t> zero_v(uint8_t(0));
for (int i = 0; i < (int)templ.features.size(); ++i)
{
Feature f = templ.features[i];
if (f.x < 0 || f.x >= size.width || f.y < 0 || f.y >= size.height)
continue;
const uchar *lm_ptr = accessLinearMemory(linear_memories, f, T, W);
int j = 0;
// *2 to avoid int8 read out of range
for(; j <= template_positions -mipp::N<int16_t>()*2; j+=mipp::N<int16_t>()){
mipp::Reg<uint8_t> src8_v((uint8_t*)lm_ptr + j);
// uchar to short, once for N bytes
mipp::Reg<int16_t> src16_v(mipp::interleavelo(src8_v, zero_v).r);
mipp::Reg<int16_t> dst_v((int16_t*)dst_ptr + j);
mipp::Reg<int16_t> res_v = src16_v + dst_v;
res_v.store((int16_t*)dst_ptr + j);
}
for(; j<template_positions; j++)
dst_ptr[j] += short(lm_ptr[j]);
}
}
similarity这里还利用了之前spread的特性，由于spread了，所以每个T x T块内就已经包含了该块的所有信息，因此只需要对下采样T之后的结果进行运算，这里也是加速。
