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

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

1、算法原理

Shi-Tomasi算法是Harris算法的改进，在Harris算法中，是根据协方差矩阵M的两个特征值的组合来判断是否角点。而在Shi-Tomasi算法中，是根据较小的特征值是否大于阈值来判断是否角点。
这个判断依据是：较小的特征值表示在该特征值方向上的方差较小，较小的方差如果都能够大于阈值，在这个方向上的变化满足是角点的判断要求。
协方差矩阵$M$如下表示：$w(x,y)$表示窗函数，通常不使用。$Ix，Iy$分别表示在x及y方向上的差分。

$$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(\lambda I-A)=0$
对于 $A$这样一个2x2矩阵来说，我们可以分解为：
$$det(\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$$
$$(\lambda-a)(\lambda-c)-bd=\lambda^2-(a+c)\lambda +ac-bd=0$$
根据求根公式，有：
$$\lambda = \frac{{(a+c)\pm \sqrt{(a+c)^2-4(ac-bd)}}}{2}=\frac{{(a+c)\pm \sqrt{(a-c)^2+4bd}}}{2}$$
带入具体的 $M$，较小的特征值为：
$$\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、先求$I_x 、I_y$；
2、再求对应的$I_x^2、I_y^2、I_xI_y$及其$\sum I_x^2$ 、$\sum I_xI_y$ 、$\sum I_y^2$
3、最后按照公式计算特征值

在计算特征值时，是调用calcMinEigenVal() 。该接口是根据以上公式

$$\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=\frac{\sum I_x^2}{2}$ $b= \sum I_xI_y$ $c=\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进行分析。