Matlab calibration board corner detection principle

Corner detection is camera calibration （ Chessboard ） A very important link in , Many people think that Matlab Calibrate the camera ratio OpenCV A more stable , One of the main reasons is its corner detection method . So , I read. Matlab Corner detection related source code （ Can be in “ The installation path /toolbox/vision/vision/+vision/+internal/+calibration/+checkerboard” Found in folder ）, Roughly summarize its principle .

Algorithm flow

No matter what corner detection method （FAST,GFTT,…）, There is a similar process ：

1. Corner response calculation ： That is to calculate the of each pixel “ score ”, To evaluate the possibility that it is a corner
2. Non maximum suppression ： In a local context , Find the corner, and the point with the largest response is recorded as corner , Other points are eliminated
3. Corner subpixelization （ Optional ）： The first 2 The corner coordinates obtained in step are integer , For application scenarios with high accuracy （ Like camera calibration ）, We need to find its sub-pixel coordinates further

Among them the first 2、3 The steps are almost the same , The non maxima is suppressed at Matlab In the corresponding find_peaks function , By calling imregionalmax Function implementation , In fact, it depends on whether the corner response of a point exceeds a specific threshold , And is the maximum value in the neighborhood , Pixels that meet these two conditions at the same time are judged as corner points . The sub-pixel corner is basically based on the neighborhood gradient voting method , Corresponding to OpenCV in cornerSubPix function .

therefore , The main factors affecting the stability of corner detection are 1 Step ： Calculation of corner response .

Corner response

Matlab Use in secondDerivCornerMetric Function to calculate the corner response of each pixel , Source code ：

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

% Low-pass filter the image
%  Gaussian filtering is used to suppress noise 
G = fspecial('gaussian', coder.const(round(sigma * 7)+1), sigma);
Ig = imfilter(I, G, 'conv');

%  One dimensional convolution kernel 
derivFilter = [-1 0 1];

% first derivatives
%  To image x and y The directions are convoluted separately , Calculate a step 
Iy = imfilter(Ig, derivFilter', 'conv');
Ix = imfilter(Ig, derivFilter, 'conv');
    
% define steerable filter constants
%  Precomputing ±45° In two directions sin/cos, As a weight for subsequent use 
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° One step in the direction 
I_45 = Ix * cosPi4 + Iy * sinPi4;
I_n45 = Ix * cosNegPi4 + Iy * sinNegPi4;

% second derivative  Two steps 
% xy Two steps in direction 
Ixy = imfilter(Ix, derivFilter', 'conv');

I_45_x = imfilter(I_45, derivFilter, 'conv');
I_45_y = imfilter(I_45, derivFilter', 'conv');
% 45° Two steps in direction 
I_45_45 = I_45_x * cosNegPi4 + I_45_y * sinNegPi4;

% suppress the outer corners
% xy Corner response in direction 
cxy = sigma^2 * abs(Ixy) - 1.5 * sigma * (abs(I_45) + abs(I_n45));
cxy(cxy < 0) = 0;
% 45° Corner response in direction 
c45 = sigma^2 * abs(I_45_45) - 1.5 * sigma * (abs(Ix) + abs(Iy));
c45(c45 < 0) = 0;

Although the theoretical explanation of this code is not found , But the idea is still relatively easy to understand . such as cxy The value of the xy The second step of the direction increases with the increase of , With 45° The direction decreases with the increase of step degree , So if a pixel is “ Ten ” Glyph corner , It's in cxy You will get a larger response value . Empathy ,c45 Able to respond “X” The response value of the corner .

The above code uses C++ be based on OpenCV Realization ：

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

Since the direction of corners in the same calibration image is consistent ,Matlab Based on cxy and c45 Go through the process of corner extraction and finding the calibration board , And calculate the energy of the calibration plate （ The search results of reaction calibration board are good or bad ）, Select the detection result with large energy as the final output result .