本文章是【opencv】goodFeaturesToTrack源码分析-1的后续,主要描述Shi-Tomasi角点检测算法原理及opencv实现。
1、算法原理
Shi-Tomasi算法是Harris算法的改进,在Harris算法中,是根据协方差矩阵M的两个特征值的组合来判断是否角点。而在Shi-Tomasi算法中,是根据较小的特征值是否大于阈值来判断是否角点。
这个判断依据是:较小的特征值表示在该特征值方向上的方差较小,较小的方差如果都能够大于阈值,在这个方向上的变化满足是角点的判断要求。
协方差矩阵
M
如下表示:
这里我们需要关心的是,如何求解矩阵M的特征值?
求解特征值,一般求解这样的一个 λ ,使得:
对于 A 这样一个2x2矩阵来说,我们可以分解为:
根据求根公式,有:
带入具体的 M ,较小的特征值为:
后续在具体的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、Iy
;
2、再求对应的
I2x、I2y、IxIy
及其
∑I2x
、
∑IxIy
、
∑I2y
3、最后按照公式计算特征值
在计算特征值时,是调用calcMinEigenVal() 。该接口是根据以上公式
从下面的代码可以看出, a=∑I2x2 b=∑IxIy c=∑I2y2
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进行分析。