Shi-Tomasi角点检测源码分析

一.Shi-Tomasi相对于Harris的改进

Shi-Tomas在角点响应函数上与Harris不同。
　1.Shi-Tomas角点响应函数
　　　　　　　 R = min(λ1，λ2)
　2.Harris角点响应函数
　　　　　　　 R = det(M) - k(Trace(M)^2)
　　其中M，是图像的张量矩阵。λ1与λ2是M的特征值。

二.Shi-Tomasi角点检测角点响应函数的优点

1.个人理解
1）在所表示的几何椭圆中，特征值越大，对应的椭圆轴越短，E(U,V)随（u,v）的改变越大。
2）当某点最小的特征值大于所设定的阈值时，矩阵M最少有两个特征值大于设定阈值。从而此点一定满足角点要求。
因此以特征值最小值作为判断点能够满足强角点 。
2.图像直观表示

三.源码分析

说明：Shi-Tomasi角点检测与我上一篇文章的Harris角点检测所分析的源码几乎一模一样，如果理解了上一篇文章，本文就非常简单了。两者唯一的不同就是最后的角点响应函数的计算公式。其余的算法流程等完全相同。

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

    CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
               ocl_cornerMinEigenValVecs(_src, _dst, blockSize, ksize, 0.0, borderType, MINEIGENVAL))

/*#ifdef HAVE_IPP
    int kerSize = (ksize < 0)?3:ksize;
    bool isolated = (borderType & BORDER_ISOLATED) != 0;
    int borderTypeNI = borderType & ~BORDER_ISOLATED;
#endif
    CV_IPP_RUN(((borderTypeNI == BORDER_REPLICATE && (!_src.isSubmatrix() || isolated)) &&
            (kerSize == 3 || kerSize == 5) && (blockSize == 3 || blockSize == 5)) && IPP_VERSION_X100 >= 800,
        ipp_cornerMinEigenVal( _src, _dst, blockSize, ksize, borderType ));
    */

    Mat src = _src.getMat();
    _dst.create( src.size(), CV_32FC1 );
    Mat dst = _dst.getMat();

    cornerEigenValsVecs( src, dst, blockSize, ksize, MINEIGENVAL, 0, borderType );//此处的模式选择为MINEIGENVAL最小特征值
}
static void
cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,//在Shi-Tomasi角点检测中，此函数调用与Harris角点检测一致，唯一的不同在最后一步的模式选择
                     int aperture_size, int op_type, double k=0.,
                     int borderType=BORDER_DEFAULT )
{
#if CV_TRY_AVX//判断AVX
    bool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
#endif
#if CV_SIMD128//判断SIMD128
    bool haveSimd = hasSIMD128();
#endif

    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 );

    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;

    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);

#if CV_TRY_AVX
        if( haveAvx )
            j = cornerEigenValsVecsLine_AVX(dxdata, dydata, cov_data, size.width);
        else
#endif // CV_TRY_AVX
            j = 0;

#if CV_SIMD128
        if( haveSimd )
        {
            for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes )
            {
                v_float32x4 v_dx = v_load(dxdata + j);
                v_float32x4 v_dy = v_load(dydata + j);

                v_float32x4 v_dst0, v_dst1, v_dst2;
                v_dst0 = v_dx * v_dx;
                v_dst1 = v_dx * v_dy;
                v_dst2 = v_dy * v_dy;

                v_store_interleave(cov_data + j * 3, v_dst0, v_dst1, v_dst2);
            }
        }
#endif // CV_SIMD128

        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;
        }
    }

    boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),
        Point(-1,-1), false, borderType );

    if( op_type == MINEIGENVAL )//此处的选择模式是最小特征值
        calcMinEigenVal( cov, eigenv );
    else if( op_type == HARRIS )
        calcHarris( cov, eigenv, k );
    else if( op_type == EIGENVALSVECS )
        calcEigenValsVecs( cov, eigenv );
}
static void calcMinEigenVal( const Mat& _cov, Mat& _dst )
{
    int i, j;
    Size size = _cov.size();
#if CV_TRY_AVX//判别AVX
    bool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
#endif
#if CV_SIMD128//判别SIMD128
    bool haveSimd = hasSIMD128();
#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);
#if CV_TRY_AVX
        if( haveAvx )
            j = calcMinEigenValLine_AVX(cov, dst, size.width);
        else
#endif // CV_TRY_AVX
            j = 0;

#if CV_SIMD128
        if( haveSimd )
        {
            v_float32x4 half = v_setall_f32(0.5f);
            for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes )
            {
                v_float32x4 v_a, v_b, v_c, v_t;
                v_load_deinterleave(cov + j*3, v_a, v_b, v_c);
                v_a *= half;
                v_c *= half;
                v_t = v_a - v_c;
                v_t = v_muladd(v_b, v_b, (v_t * v_t));
                v_store(dst + j, (v_a + v_c) - v_sqrt(v_t));
            }
        }
#endif // CV_SIMD128
		//直接进行串行处理

        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));//2阶矩阵两个特征值分别是
																		//λ1 = (float)((a + c) + std::sqrt((a - c)*(a - c) + b*b))
        }																//λ2 = (float)((a + c) - std::sqrt((a - c)*(a - c) + b*b))这个是最小的
    }
}