shi-Tomasi score

发表于 | 分类于

在视觉里程计部分，在对图像帧提取Fast特征点后，计算shi-Tomasi角点检测，根据阈值选择更好的角点。

那什么是shi-Tomasi角点检测？shi-Tomasi分数又该如何计算呢？shi-Tomasi是harris角点的改进版。我们先来看看harris角点检测。

harris 角点检测

harris的想法是，在image里，如果一个小的patch（或者叫窗口）进行小范围的移动时，窗口内的像素灰度值能产生明显的变化，则认为该窗口中存在角点。转换成数学表达式为：

  • $E$是原始窗口和移动后的窗口的差分
  • $u$是在$x$方向上窗口的平移
  • $v$是在$y$方向上窗口的平移
  • $w(x,y)$是窗口在点$(x,y)$位置上的权重
  • $I$是图像在点$(x,y)$处的灰度值
  • $I(x+u,y+v)$指窗口移动后的灰度值
  • $I(x,y)$指原始窗口的灰度值

现在若要找到一个harris角点，则需要寻找一个窗口，使得其能产生一个大的$E$。也就是要最大化下面的式子(这里先把权重放一边，后面再补上)：

对上式进行一阶泰勒展开得：

在这里把刚才我们丢掉的权重补上，令

所以，

$M$的特征值$\lambda_1,\lambda_2$有助于判断一个窗口是否含有角点：

  • 若两个特征值都小，则在各个方向上都没有变化，说明窗口处在平滑区域
  • 若其中一个特征值大，另一个小，说明窗口在边缘上，沿边缘方向移动没有变化
  • 若两个特征值都大，则窗口中存在角点

至此，我们就可以计算harris的评分分数$R$，若分数$R$大于某一阈值，则认为该窗口中存在角点：

shi-Tomasi角点检测

shi-Tomasi角点检测是在harris上稍作修改，修改了评分计算：

在视觉里程计中，该部分的实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
float FastDetector::shiTomasiScore(const cv::Mat& img, int u, int v)
	{
		assert(img.type() == CV_8UC1);

		float dXX = 0.0;
		float dYY = 0.0;
		float dXY = 0.0;
		const int halfbox_size = 4;
		const int box_size = 2 * halfbox_size;
		const int box_area = box_size*box_size;
		const int x_min = u - halfbox_size;
		const int x_max = u + halfbox_size;
		const int y_min = v - halfbox_size;
		const int y_max = v + halfbox_size;

		if (x_min < 1 || x_max >= img.cols - 1 || y_min < 1 || y_max >= img.rows - 1)
			return 0.0; // 面片太靠近边界，返回0

		const int stride = img.step.p[0];//一行元素的个数
		for (int y = y_min; y < y_max; ++y)//垂直平移,计算每次平移后对角线上的灰度值差
		{
			const uint8_t* ptr_left = img.data + stride*y + x_min - 1;
			const uint8_t* ptr_right = img.data + stride*y + x_min + 1;
			const uint8_t* ptr_top = img.data + stride*(y - 1) + x_min;
			const uint8_t* ptr_bottom = img.data + stride*(y + 1) + x_min;
			for (int x = 0; x < box_size; ++x, ++ptr_left, ++ptr_right, ++ptr_top, ++ptr_bottom)
			{
				float dx = *ptr_right - *ptr_left;//计算I在x处的导数
				float dy = *ptr_bottom - *ptr_top;//计算I在y处的导数
				dXX += dx*dx;
				dYY += dy*dy;
				dXY += dx*dy;
			}
		}

		// 返回小的特征值
		dXX = dXX / (2.0 * box_area);//应该是乘了w(x,y)，权重
		dYY = dYY / (2.0 * box_area);
		dXY = dXY / (2.0 * box_area);
    //解一元二次方程，求最小特征值
		return 0.5 * (dXX + dYY - sqrt((dXX + dYY) * (dXX + dYY) - 4 * (dXX * dYY - dXY * dXY)));
	}