【opencv】goodFeaturesToTrack源码分析-2-Shi-Tomasi角点检测

本文章是【opencv】goodFeaturesToTrack源码分析-1的后续，主要描述Shi-Tomasi角点检测算法原理及opencv实现。

1、算法原理

Shi-Tomasi算法是Harris算法的改进，在Harris算法中，是根据协方差矩阵M的两个特征值的组合来判断是否角点。而在Shi-Tomasi算法中，是根据较小的特征值是否大于阈值来判断是否角点。

这个判断依据是：较小的特征值表示在该特征值方向上的方差较小，较小的方差如果都能够大于阈值，在这个方向上的变化满足是角点的判断要求。

协方差矩阵MM如下表示：w(x,y)w(x,y)表示窗函数，通常不使用。Ix，IyIx，Iy分别表示在x及y方向上的差分。

M=∑x,yw(x,y)[I2xIxIyIxIyI2y]M= \sum_{x,y}w(x,y) \left[ \begin{matrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{matrix} \right]

这里我们需要关心的是，如何求解矩阵M的特征值？

求解特征值，一般求解这样的一个λ\lambda，使得:

det(λI−A)=0det(\lambda I-A)=0

对于AA这样一个2x2矩阵来说，我们可以分解为：

det(λ[1001]−A)=det([λ00λ]−[adbc])=0det(\lambda \left[ \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}
\right]-A
)= det( \left[ \begin{matrix} \lambda& 0 \\ 0 & \lambda \end{matrix}
\right]-
\left[
\begin{matrix}
a &b \\
d& c
\end{matrix}
\right]
)=0

(λ−a)(λ−c)−bd=λ2−(a+c)λ+ac−bd=0(\lambda-a)(\lambda-c)-bd=\lambda^2-(a+c)\lambda +ac-bd=0

根据求根公式，有：

λ=(a+c)±(a+c)2−4(ac−bd)−−−−−−−−−−−−−−−−−√2=(a+c)±(a−c)2+4bd−−−−−−−−−−−√2\lambda = \frac{{(a+c)\pm \sqrt{(a+c)^2-4(ac-bd)}}}{2}=\frac{{(a+c)\pm \sqrt{(a-c)^2+4bd}}}{2}

带入具体的MM，较小的特征值为：

λ=(∑I2x+∑I2y)−(∑I2x−∑I2y)2+4(∑IxIy)2−−−−−−−−−−−−−−−−−−−−−−−−√2\lambda =\frac{{(\sum I_x^2+\sum I_y^2)- \sqrt{(\sum I_x^2-\sum I_y^2)^2+4(\sum I_xI_y)^2}}}{2}

后续在具体的opencv源码分析中就会看到该公式的应用。

2、算法实现

该算法是作为计算最小特征值被goodFeaturesToTrack调用的，调用的位置为：

void cv::goodFeaturesToTrack( )
{
...
// 计算最小特征值
cornerMinEigenVal( image, eig, blockSize, 3 );
...
}

在：..\sources\modules\imgproc\src\Featureselect.cpp中。

可以看到该算法的接口是

void cv::cornerMinEigenVal( InputArray _src, OutputArray _dst, int blockSize, int ksize, int borderType ){

cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,
int aperture_size, int op_type, double k=0.,
int borderType=BORDER_DEFAULT )
}

具体实现在：…\sources\modules\imgproc\src\Corner.cpp

static void
cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,
int aperture_size, int op_type, double k=0.,
int borderType=BORDER_DEFAULT )
{
//计算缩放参数
int depth = src.depth();
double scale = (double)(1 << ((aperture_size > 0 ? aperture_size : 3) - 1)) * block_size;
if( aperture_size < 0 )
scale *= 2.0;
if( depth == CV_8U )
scale *= 255.0;
scale = 1.0/scale;

CV_Assert( src.type() == CV_8UC1 || src.type() == CV_32FC1 );

//sobel或者scharr计算分别在x方向及y方向的梯度，即Ix Iy
Mat Dx, Dy;
if( aperture_size > 0 )
{
Sobel( src, Dx, CV_32F, 1, 0, aperture_size, scale, 0, borderType );
Sobel( src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType );
}
else
{
Scharr( src, Dx, CV_32F, 1, 0, scale, 0, borderType );
Scharr( src, Dy, CV_32F, 0, 1, scale, 0, borderType );
}

Size size = src.size();
Mat cov( size, CV_32FC3 );
int i, j;
// 计算Ix^2 Iy^2 IxIy
for( i = 0; i < size.height; i++ )
{
float* cov_data = cov.ptr<float>(i);
const float* dxdata = Dx.ptr<float>(i);
const float* dydata = Dy.ptr<float>(i);
j = 0;

#if CV_NEON
for( ; j <= size.width - 4; j += 4 )
{
float32x4_t v_dx = vld1q_f32(dxdata + j);
float32x4_t v_dy = vld1q_f32(dydata + j);

float32x4x3_t v_dst;
v_dst.val[0] = vmulq_f32(v_dx, v_dx);
v_dst.val[1] = vmulq_f32(v_dx, v_dy);
v_dst.val[2] = vmulq_f32(v_dy, v_dy);

vst3q_f32(cov_data + j * 3, v_dst);
}
for( ; j < size.width; j++ )
{
float dx = dxdata[j];
float dy = dydata[j];

cov_data[j*3] = dx*dx;
cov_data[j*3+1] = dx*dy;
cov_data[j*3+2] = dy*dy;
}
}

//求 ∑Ix^2 ∑Iy^2 ∑IxIy
boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),
Point(-1,-1), false, borderType );

// 调用 calcMinEigenVal 计算特征值
if( op_type == MINEIGENVAL )
calcMinEigenVal( cov, eigenv );
else if( op_type == HARRIS )
calcHarris( cov, eigenv, k );
else if( op_type == EIGENVALSVECS )
calcEigenValsVecs( cov, eigenv );
}

在上述代码中，我们可以看到，代码实现是按照算法原理来完成的:

1、先求Ix、IyI_x 、I_y；

2、再求对应的I2x、I2y、IxIyI_x^2、I_y^2、I_xI_y及其∑I2x \sum I_x^2 、∑IxIy \sum I_xI_y 、∑I2y \sum I_y^2

3、最后按照公式计算特征值

在计算特征值时，是调用calcMinEigenVal() 。该接口是根据以上公式λ=(∑I2x+∑I2y)−(∑I2x−∑I2y)2+4(∑IxIy)2−−−−−−−−−−−−−−−−−−−−−−−−√2=(∑I2x2+∑I2y2)−(∑I2x2−∑I2y2)2+(∑IxIy)2−−−−−−−−−−−−−−−−−−−−−−−−−√\lambda =\frac{{(\sum I_x^2+\sum I_y^2)- \sqrt{(\sum I_x^2-\sum I_y^2)^2+4(\sum I_xI_y)^2}}}{2}= {(\frac{\sum I_x^2}{2}+\frac{\sum I_y^2}{2})- \sqrt{(\frac{\sum I_x^2}{2}-\frac{\sum I_y^2}{2})^2+(\sum I_xI_y)^2}} 计算特征值。

从下面的代码可以看出，a=∑I2x2a=\frac{\sum I_x^2}{2} b=∑IxIyb= \sum I_xI_y c=∑I2y2c=\frac{\sum I_y^2}{2}

static void calcMinEigenVal( const Mat& _cov, Mat& _dst )
{
int i, j;
Size size = _cov.size();
#if CV_SSE
volatile bool simd = checkHardwareSupport(CV_CPU_SSE);
#endif

if( _cov.isContinuous() && _dst.isContinuous() )
{
size.width *= size.height;
size.height = 1;
}

for( i = 0; i < size.height; i++ )
{
const float* cov = _cov.ptr<float>(i);
float* dst = _dst.ptr<float>(i);
j = 0;

#if CV_NEON
float32x4_t v_half = vdupq_n_f32(0.5f);
for( ; j <= size.width - 4; j += 4 )
{
float32x4x3_t v_src = vld3q_f32(cov + j * 3);
float32x4_t v_a = vmulq_f32(v_src.val[0], v_half);
float32x4_t v_b = v_src.val[1];
float32x4_t v_c = vmulq_f32(v_src.val[2], v_half);

float32x4_t v_t = vsubq_f32(v_a, v_c);
v_t = vmlaq_f32(vmulq_f32(v_t, v_t), v_b, v_b);
vst1q_f32(dst + j, vsubq_f32(vaddq_f32(v_a, v_c), cv_vsqrtq_f32(v_t)));
}
#endif
for( ; j < size.width; j++ )
{
float a = cov[j*3]*0.5f;
float b = cov[j*3+1];
float c = cov[j*3+2]*0.5f;
dst[j] = (float)((a + c) - std::sqrt((a - c)*(a - c) + b*b));
}
}
}

值得注意的是opencv中也提供NEON优化的代码，根据是否定义宏CV_NEON来判断是否使用NEON优化。对于学习NEON优化很有启发。

后续将对这里的sobel滤波及其boxfilter进行分析。