cornerSubPixel亚像素角点检测源码分析（最小二乘法）

OpenCV的中有cornerSubPixel（）这个API函式用来针对初始的整数角点坐标进行亚像素精度的优化，该函式原型如下：

C++: void cornerSubPix(InputArray image, InputOutputArray corners, Size winSize, Size zeroZone, TermCriteria criteria);

image为输入的单通道影象; corners为提取的初始整数角点;winsize为求取亚像素角点的视窗大小，比如设定尺寸（11,11），需要注意的是11为半径，则视窗大小为23x23; zeroZone是设定的“零区域”，在搜寻视窗内，设定的“零区域”内的值不会被累加，权重值为0。如果设定为大小（-1，-1），则表示没有这样的区域; critteria是条件阈值，包括迭代次数阈值和误差精度阈值，一旦其中一项条件满足设定的阈值，则停止迭代，获得亚画素角点。

这个API通过下面示例的语句调用：

cv::cornerSubPix(grayImg, pts, cv::Size(11, 11), cv::Size(-1, -1), cv::TermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 30, 0.1));

首先看标准包含的两个条件阈值在源码中是怎么设定的。如下所示，默认最大迭代次数为100次，误差精度为EPS * EPS，默认0.1×0.1。

const  int MAX_ITERS = 100 ;
int win_w = win.width * 2 + 1，win_h = win.height * 2 + 1 ;
int i，j，k;
int max_iters =（criteria.type＆CV_TERMCRIT_ITER）？MIN（MAX（criteria.maxCount，1 ），MAX_ITERS）：MAX_ITERS;
double eps =（criteria.type＆CV_TERMCRIT_EPS）？MAX（criteria.epsilon，0.）：0 ;
eps * = eps; // 在比较操作中使用平方误差

然后是高斯权重的计算，如下所示，视窗中心附近权重高，越往视窗边界权重越小如果设定的有“零区域”，则权重值设定为0。计算出的权重分布如下图：
　　Mat maskm(win_h, win_w, CV_32F), subpix_buf(win_h+2, win_w+2, CV_32F);
float* mask = maskm.ptr();

for( i = 0; i < win_h; i++ )
{
    float y = (float)(i - win.height)/win.height;
    float vy = std::exp(-y*y);
    for( j = 0; j < win_w; j++ )
    {
        float x = (float)(j - win.width)/win.width;
        mask[i * win_w + j] = (float)(vy*std::exp(-x*x));
    }
}

// make zero_zone
if( zeroZone.width >= 0 && zeroZone.height >= 0 &&
    zeroZone.width * 2 + 1 < win_w && zeroZone.height * 2 + 1 < win_h )
{
    for( i = win.height - zeroZone.height; i <= win.height + zeroZone.height; i++ )
    {
        for( j = win.width - zeroZone.width; j <= win.width + zeroZone.width; j++ )
        {
            mask[i * win_w + j] = 0;
        }
    }
}
![在这里插入图片描述](https://img-blog.csdnimg.cn/20181211160242746.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMwMzM5NTk1,size_16,color_FFFFFF,t_70)
接下来就是针对每个初始角点，按照上述公式，逐个进行迭代求取亚像素角点，源码如下。

// do optimization loop for all the points
for( int pt_i = 0; pt_i < count; pt_i++ )
{
Point2f cT = corners[pt_i], cI = cT;
int iter = 0;
double err = 0;

    do
    {
        Point2f cI2;
        double a = 0, b = 0, c = 0, bb1 = 0, bb2 = 0;

        getRectSubPix(src, Size(win_w+2, win_h+2), cI, subpix_buf, subpix_buf.type());
        const float* subpix = &subpix_buf.at<float>(1,1);

        // process gradient
        for( i = 0, k = 0; i < win_h; i++, subpix += win_w + 2 )
        {
            double py = i - win.height;

            for( j = 0; j < win_w; j++, k++ )
            {
                double m = mask[k];
                double tgx = subpix[j+1] - subpix[j-1];
                double tgy = subpix[j+win_w+2] - subpix[j-win_w-2];
                double gxx = tgx * tgx * m;
                double gxy = tgx * tgy * m;
                double gyy = tgy * tgy * m;
                double px = j - win.width;

                a += gxx;
                b += gxy;
                c += gyy;

                bb1 += gxx * px + gxy * py;
                bb2 += gxy * px + gyy * py;
            }
        }

        double det=a*c-b*b;
        if( fabs( det ) <= DBL_EPSILON*DBL_EPSILON )
            break;

        // 2x2 matrix inversion
        double scale=1.0/det;
        cI2.x = (float)(cI.x + c*scale*bb1 - b*scale*bb2);
        cI2.y = (float)(cI.y - b*scale*bb1 + a*scale*bb2);
        err = (cI2.x - cI.x) * (cI2.x - cI.x) + (cI2.y - cI.y) * (cI2.y - cI.y);
        cI = cI2;
        if( cI.x < 0 || cI.x >= src.cols || cI.y < 0 || cI.y >= src.rows )
            break;
    }
    while( ++iter < max_iters && err > eps );

    // if new point is too far from initial, it means poor convergence.
    // leave initial point as the result
    if( fabs( cI.x - cT.x ) > win.width || fabs( cI.y - cT.y ) > win.height )
        cI = cT;

    corners[pt_i] = cI;
}

①程式码中CI2为本次迭代获取的亚像素位置，CI为上次迭代获取的亚像素角点位置，CT是初始的整数角点位置。

②每次迭代结束计算CI与CI2之间的欧式距离err，如果两者之间的欧式距离err小于设定的阈值，或者迭代次数达到设定的阈值，则停止迭代。

③停止迭代后，需要再次判断最终的亚像素角点位置和初始整数角点之间的差异，如果差值大于设定视窗尺寸的一半，则说明最小二乘计算中收敛性不好，丢弃计算得到的亚像素角点，仍然使用初始的整数角点。
　 ps:要理解代码请看：
　
代入推导得到的公式即可。
原文请看https://tw.saowen.com/a/6ffb4ed611d5582399844d286fbe93f58a745ca7731b71bda7de08d613124b87