OpenCV Shi-Tomasi角点解析笔记

会飞的鱼chelmx

已于 2023-03-23 00:40:19 修改

阅读量193

点赞数

分类专栏： SLAM 文章标签： opencv 计算机视觉人工智能

于 2023-02-05 18:01:43 首次发布

本文链接：https://blog.csdn.net/m0_37142194/article/details/128892245

版权

SLAM 专栏收录该内容

8 篇文章 1 订阅

订阅专栏

文章目录

1. Shi-Tomasi角点数学原理

1.1 角点的定义

首先通过固定大小的窗口观察图像，然后使窗口按任意方向滑动再次观察图像，比较滑动前后两种情况窗口中的像素灰度变化程度，如果在任意方向上滑动，像素灰度都有较大的变化，那么认为该窗口中存在角点（corner）。

下图给出了三种常见情况：

$\text{flat}$ ：该窗口在任何方向上滑动，灰度变化程度不大。
$\text{edge}$ ：该窗口在某个特定方向上滑动，灰度变化程度不大。
$\text{corner}$ ：该窗口在任何方向上滑动，灰度变化程度明显。

在这里插入图片描述

1.2 角点的数学推导

当窗口按 $[u, v]$ 方向进行移动时，滑动前后像素灰度变化程度可用下式表示：
$E(u,v)=\sum_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)]^2$
其中 $w(\cdot)$ 表示窗口函数， $I(\cdot)$ 表示图像函数。

将 $I (x + u, y + v)$ 利用一阶泰勒展开有：
$I(x+u,y+v)\approx I(x,y)+I_xu+I_yv$
其中 $I_x$ 和 $I_y$ 表示沿图像沿 $x$ 和 $y$ 方向的一阶导数，可用 $\text{Sobel}$ 算子计算。

代入泰勒展开则有：
$\begin{aligned} E(u,v)&\approx\sum_{x,y}w(x,y)[I(x,y)+I_xu+I_yv-I(x,y)]^2\\ &\approx\sum_{x,y}w(x,y)[I_x^2u^2+2I_xI_yuv+I_y^2v^2]\\ &\approx\sum_{x,y}w(x,y)\begin{bmatrix}u&v\end{bmatrix}\begin{bmatrix}I_x^2&I_xI_y\\I_xI_y&I_y^2\end{bmatrix}\begin{bmatrix}u\\v\end{bmatrix}\\ &\approx\begin{bmatrix}u&v\end{bmatrix}\sum_{x,y}w(x,y)\begin{bmatrix}I_x^2&I_xI_y\\I_xI_y&I_y^2\end{bmatrix}\begin{bmatrix}u\\v\end{bmatrix}\\ &\approx\begin{bmatrix}u&v\end{bmatrix}\begin{bmatrix}\sum_{x,y}w(x,y)I_x^2&\sum_{x,y}w(x,y)I_xI_y\\\sum_{x,y}w(x,y)I_xI_y&\sum_{x,y}w(x,y)I_y^2\end{bmatrix}\begin{bmatrix}u\\v\end{bmatrix}\\ &\approx\begin{bmatrix}u&v\end{bmatrix}M\begin{bmatrix}u\\v\end{bmatrix} \end{aligned}$
上式右侧的表达式为一个二项式函数，其本质是一个椭圆函数：
$Au^2+2Buv+Cv^2=1$
其中 $A$ 表示 $m_{11}$ ， $B$ 表示 $m_{12}=m_{21}$ ， $C$ 表示 $m_{22}$ 。椭圆的扁率和尺寸由 $M$ 的特征值 $\lambda_1$ 、 $\lambda_2$ 决定，椭圆的方向由 $M$ 的特征矢量决定。

特征根与椭圆关系

根据上述分析， $M$ 特征值大小比值可分为三种情况：

一个特征值大，另一个特征值小， $\lambda_1\gg\lambda_2$ 或 $\lambda_2\gg\lambda_1$ 。说明窗口在某个特定方向上滑动，灰度变化程度不大，在垂直方向上滑动，灰度变化程度明显，可以判断窗口存在** $\text{edge}$ **。
两个特征值都小，且近似相等。说明窗口在任何方向上滑动，灰度变化程度不大，可以判断窗口为 $\text{flat}$ 区域。
两个特征值都大，且近似相等。说明窗口在任何方向上滑动，灰度变化程度明显，可以判断窗口存在 $\text{corner}$ 。

特征根三种情况

1.3 Shi-Tomasi角点的判定

对于 $\text{Shi-Tomasi}$ 角点而言，其判定条件不要求两特征值近似相等，仅要求两特征值较小的一个大于特定的阈值。

根据 $2\times2$ 矩阵的特征值推导公式可知，较小的特征值为：
$\begin{aligned} \lambda&=\frac{(A+C)-\sqrt{(A-C)^2+4B^2}}{2}\\ &=\frac{(\sum_{x,y}w(x,y)I_x^2+\sum_{x,y}w(x,y)I_y^2)-\sqrt{(\sum_{x,y}w(x,y)I_x^2-\sum_{x,y}w(x,y)I_y^2)+4(\sum_{x,y}w(x,y)I_xI_y)^2}}{2} \end{aligned}$

2. Shi-Tomasi角点代码解析

在 $\text{OpenCV}$ 中，关于 $\text{Shi-Tomasi}$ 算子的核心代码位于 $\text{corner.cpp}$ 文件中：

void cv::cornerMinEigenVal( InputArray _src, OutputArray _dst, int blockSize, int ksize, int borderType )
                                                                // _src：输入图像
                                                                // _dst：输出角点
                                                                // blockSize：观察窗尺寸
                                                                // ksize：一阶差分算子阶数
                                                                // borderType：边缘处理类型
{
    CV_INSTRUMENT_REGION();

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

/*#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 );
}

static void
cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,
                     int aperture_size, int op_type, double k=0.,
                     int borderType=BORDER_DEFAULT )
                                                                // src：输入图像
                                                                // eigenv：输出特征根
                                                                // block_size：观察窗尺寸
                                                                // aperture_size：一阶差分算子阶数
                                                                // op_type：操作类型
                                                                // k = 0.：Harris 角点系数
                                                                // borderType = BORDER_DEFAULT：边缘处理类型
                                                                // ktype：算子数据类型
{
#if CV_TRY_AVX
    bool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
#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 );
                                                                // 生成x方向Sobel滤波结果（一阶差分）
        Sobel( src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType );
                                                                // 生成y方向Sobel滤波结果（一阶差分）
    }
    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
        {
            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;                              // 卷积结果通道1保存x方向一阶差分平方结果
            cov_data[j*3+1] = dx*dy;                            // 卷积结果通道2保存xy方向一阶差分相乘结果
            cov_data[j*3+2] = dy*dy;                            // 卷积结果通道3保存y方向一阶差分平方结果
        }
    }

    boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),
        Point(-1,-1), false, borderType );                      // 盒式滤波
                                                                // cov：输入图像
                                                                // cov：输出图像
                                                                // cov.depth()：输出图像矩阵类型
                                                                // Size(block_size, block_size)：核大小
                                                                // Point(-1, -1)：内核中的锚点位置，默认值 (-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 )
                                                                // _cov：输入特征图，通道1保存x方向一阶差分平方结果，通道2保存xy方向一阶差分相乘结果，通道3保存y方向一阶差分平方结果
                                                                // _dst：输出特征根
{
    int i, j;
    Size size = _cov.size();
#if CV_TRY_AVX
    bool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
#endif

    if( _cov.isContinuous() && _dst.isContinuous() )            // 如果_cov及_dst连续存储，则调整size大小以加速以下循环代码
    {
        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
        {
            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));
                                                                // 求较小特征值
        }
    }
}