Matlab标定板角点检测原理

角点检测是相机标定（棋盘格）中很重要的一个环节，很多人认为用 Matlab 标定相机比 OpenCV 更稳定，主要的原因之一便在于其角点检测方法。为此，我阅读了 Matlab 中角点检测相关源码（可以在“安装路径/toolbox/vision/vision/+vision/+internal/+calibration/+checkerboard”文件夹中找到），大致总结一下其原理。

算法流程

不管是什么角点检测方法（FAST，GFTT，…），都有一个类似的流程：

1. 角点响应计算：也就是计算每个像素点的“得分”，来评价其是一个角点的可能性
2. 非极大值抑制：在一个局部范围内，找到角点响应最大的点记为角点，其它的点则被剔除
3. 角点亚像素化（可选）：第2步得到的角点坐标是整数的，对于精度较高的应用场景（比如相机标定），还需要进一步查找其亚像素坐标

其中第2、3步都差不多，非极大值抑制在 Matlab 中对应 find_peaks 函数，通过调用 imregionalmax 函数实现，其实也就是看一个点的角点响应是否超过特定阈值，并且是邻域范围内的最大值，同时满足这两个条件的像素点被判为角点。角点亚像素化采用的基本上都是邻域梯度投票法，对应于 OpenCV 中 cornerSubPix 函数。

因此，影响角点检测稳定性的主要是第1步：角点响应的计算。

角点响应

Matlab 中使用 secondDerivCornerMetric 函数计算每个像素的角点响应，源码：

function [cxy, c45, Ix, Iy, Ixy, I_45_45] = secondDerivCornerMetric(I, sigma)
%#codegen

% Low-pass filter the image
% 使用高斯滤波抑制噪声
G = fspecial('gaussian', coder.const(round(sigma * 7)+1), sigma);
Ig = imfilter(I, G, 'conv');

% 一维卷积核
derivFilter = [-1 0 1];

% first derivatives
% 对图像x和y方向分别进行卷积，计算一阶梯度
Iy = imfilter(Ig, derivFilter', 'conv');
Ix = imfilter(Ig, derivFilter, 'conv');
    
% define steerable filter constants
% 预先计算±45°两个方向的sin/cos，作为后续使用的权重
cosPi4 = coder.const(cast(cos(pi/4), 'like', I));
cosNegPi4 = coder.const(cast(cos(-pi/4), 'like', I));
sinPi4 = coder.const(cast(sin(pi/4), 'like', I));
sinNegPi4 = coder.const(cast(sin(-pi/4), 'like', I));

% first derivative at 45 degrees
% 45°方向的一阶梯度
I_45 = Ix * cosPi4 + Iy * sinPi4;
I_n45 = Ix * cosNegPi4 + Iy * sinNegPi4;

% second derivative 二阶梯度
% xy方向的二阶梯度
Ixy = imfilter(Ix, derivFilter', 'conv');

I_45_x = imfilter(I_45, derivFilter, 'conv');
I_45_y = imfilter(I_45, derivFilter', 'conv');
% 45°方向的二阶梯度
I_45_45 = I_45_x * cosNegPi4 + I_45_y * sinNegPi4;

% suppress the outer corners
% xy方向的角点响应
cxy = sigma^2 * abs(Ixy) - 1.5 * sigma * (abs(I_45) + abs(I_n45));
cxy(cxy < 0) = 0;
% 45°方向的角点响应
c45 = sigma^2 * abs(I_45_45) - 1.5 * sigma * (abs(Ix) + abs(Iy));
c45(c45 < 0) = 0;

虽然没有找到这段代码的理论说明，但其思路还是比较容易理解的。比如 cxy 的值随着 xy 方向的二阶梯度的增大而增大，随着45°方向一阶梯度的增大而减小，所以如果一个像素点是“十”字形角点，它在 cxy 中就会获得较大的响应值。同理，c45 能够反应“X”形角点的响应值。

上述代码用 C++ 基于 OpenCV 实现：

void secondDerivCornerMetric(const cv::Mat& I, cv::Mat& cxy, cv::Mat& c45,
    cv::Mat& Ix, cv::Mat& Iy, cv::Mat& Ixy, cv::Mat& I_45_45)
{
    using Scalar = double;
    constexpr int DEPTH = cv::DataType<Scalar>::depth;
    auto sigma = Scalar(2.0);
    // Low-pass filter the image
    cv::Mat Ig;
    {
        cv::Mat _I;
        I.convertTo(_I, DEPTH);
        int d = (cvRound(sigma * 7) + 1) | 1;
        GaussianBlur(_I, Ig, cv::Size(d, d), sigma);
    }

    // first derivatives
    const cv::Matx13d derivFilter(1, 0, -1);
    filter2D(Ig, Ix, DEPTH, derivFilter);
    filter2D(Ig, Iy, DEPTH, derivFilter.t());

    // define steerable filter constants
    constexpr float cosPi4 = 0.707106781186548; // cos(pi/4)
    constexpr float cosNegPi4 = cosPi4;
    constexpr float sinPi4 = cosPi4;
    constexpr float sinNegPi4 = -cosPi4;

    // first derivative at 45 degrees
    cv::Mat I_45 = Ix * cosPi4 + Iy * sinPi4;
    cv::Mat I_n45 = Ix * cosNegPi4 + Iy * sinNegPi4;

    // second derivative
    filter2D(Ix, Ixy, DEPTH, derivFilter.t());
    cv::Mat I_45_x, I_45_y;
    filter2D(I_45, I_45_x, DEPTH, derivFilter);
    filter2D(I_45, I_45_y, DEPTH, derivFilter.t());
    I_45_45 = I_45_x * cosNegPi4 + I_45_y * sinNegPi4;

    // suppress the outer corners
    auto removeNegative = [](auto& pix, const int*) { if( pix < 0 ) pix = 0; };
    cxy = sigma * sigma * abs(Ixy) - 1.5 * sigma * (abs(I_45) + abs(I_n45));
    cxy.forEach<Scalar>(removeNegative);
    c45 = sigma * sigma * abs(I_45_45) - 1.5 * sigma * (abs(Ix) + abs(Iy));
    c45.forEach<Scalar>(removeNegative);
}

由于同一张标定图像中的角点方向是一致的，Matlab 中分别基于 cxy 和 c45 走一遍角点提取和查找标定板的流程，并计算标定板能量（反应标定板查找结果的好坏），选择能量较大的检测结果作为最终输出结果。