基于生长的棋盘格角点检测方法--（2）代码详解（上）

上一篇介绍了基于生长的棋盘格角点检测方法的大概原理，详见：基于生长的棋盘格角点检测方法–（1）原理介绍
本文进一步从代码解读角度出发，更深入地理解工程中是如何实现的。
本文中用到的代码可以从以下链接下载
http://www.cvlibs.net/software/libcbdetect/
这里我把代码中主要的函数提取出来作为算法骨架，这样比较好和论文对应，可以帮助读者在茫茫代码中抓住重点。
代码框架结构如下，包括了主要的函数。其中缩进表示包含从属关系。

代码框架

I = imread('image.jpg');
corners = findCorners(I,0.01,1);
    function template = createCorrelationPatch(angle_1,angle_2,radius)
    corners.p = nonMaximumSuppression(img_corners,3,0.025,5);
    corners = refineCorners(img_du,img_dv,img_angle,img_weight,corners,10);
        [v1,v2] = edgeOrientations(img_angle_sub,img_weight_sub);
    corners = scoreCorners(img,img_angle,img_weight,corners,radius);
chessboards = chessboardsFromCorners(corners);
    chessboard = initChessboard(corners,i);
    energy = chessboardEnergy(chessboard,corners)
    proposal{j} = growChessboard(chessboard,corners,j);
        pred = predictCorners(p1,p2,p3)
        idx = assignClosestCorners(cand,pred);
plotChessboards(chessboards,corners);

本篇先介绍第一个重要函数：

findCorners

该函数的目的是从一幅包含棋盘（可以是多个）自然图像中找到棋盘中每个角点的位置。首先利用自定义的模板来突出角点，效果类似于显著性检测。然后用非极大值抑制算法来获得极大值候选点。然后对这些候选点进行亚像素级精细化（refinement），最后根据一定的规则对每个角点进行评分，最后得到较为纯净的角点。

createCorrelationPatch

首先就是创建模板（即论文中所说的prototypes），实际使用中考虑了图像中棋盘尺度的不同尺寸，所以分别创建了3个不同尺度的模板原型，实现代码如下

template_props = [0 pi/2 radius(1); pi/4 -pi/4 radius(1); 0 pi/2 radius(2); pi/4 -pi/4 radius(2); 0 pi/2 radius(3); pi/4 -pi/4 radius(3)];
function template = createCorrelationPatch(angle_1,angle_2,radius)

得到的3个不同尺度下的模板原型如下：

尺度1


尺度2


尺度3


以上三种尺度基本可以处理实际图像中所有尺度的棋盘。然后，如下代码

  img_corners_a1 = conv2(img,template.a1,'same');
  img_corners_a2 = conv2(img,template.a2,'same');
  img_corners_b1 = conv2(img,template.b1,'same');
  img_corners_b2 = conv2(img,template.b2,'same');

分别对应 ，根据论文中公式（1）可以计算出 （即case 1: a=white, b=black）
公式（1）如下，具体解释见上一篇博客。

之所以分case1，case2（对应 ）是由于棋盘格根据黑白分布顺序不同有两种，分别是
b w w b
w b b w
其中,b=black, w=white
按照上面的公式（1）即可计算每个像素的Corner likelihood。

nonMaximumSuppression

非极大值抑制（NMS）采用的窗口范围是(n+1)*(n+1), n=3。经过非极大值抑制后得到的候选点位置：

Img_weight是x,y方向梯度的2范数，后面用来作为方向直方图的加权。其结果如下图

refineCorner

对于每个Corner，以该像素点坐标为中心，取21*21窗口，代码如下

img_angle_sub  = img_angle(max(cv-r,1):min(cv+r,height),max(cu-r,1):min(cu+r,width));
img_weight_sub = img_weight(max(cv-r,1):min(cv+r,height),max(cu-r,1):min(cu+r,width));

img_angle，img_weight分别是像素点的梯度方向角和幅值。

edgeOrientations

函数 edgeOrientations 首先将窗口内所有的梯度方向映射到一个32bin的直方图里，用梯度幅值作为加权值，代码如下：

    bin = max(min(floor(vec_angle(i)/(pi/bin_num)),bin_num-1),0)+1;
    angle_hist(bin) = angle_hist(bin)+vec_weight(i);

然后用meanshift方法来寻找该直方图的局部极大值。首先要先对直方图做一个高斯平滑。

findModesMeanShift函数目的就是寻找直方图中两个最大的峰值位置，对应最大的两个梯度。正常角点的直方图应该是下图左，有两个幅值相当的峰值。若两个峰值差别太大（下图右），则认为该点不是角点。

我们把这两个最大峰值对应的位置称为该直方图的modes。然后计算这两个modes直接对应的梯度角的差值，这个差值要大于一定的阈值才认为该角点有效。
然后就是corner location refinement

        w  = [u v]-[cu cv];
        d1 = norm(w-w*v1'*v1);
        d2 = norm(w-w*v2'*v

其中w*v1’*v1表示向量 w在单位向量v1方向的投影向量。

refineCorners函数结束后，滤掉不合格的Corner。如下图，绿色表示refine后剩余的Corner，红色表示滤掉的伪Corner

ScoreCorners

取出角点的radius* radius邻域（仍然分3个不同尺度），根据当前角点的两个主方向向量创建模板。
score_intensity的计算参考论文（见最后的参考文献）中公式（1），和前面计算Corner likelihood的方法一样。
最后的score是梯度得分和likelihood得分的乘积。

template=reateCorrelationPatch(atan2(v1(2),v1(1)),atan2(v2(2),v2(1)),c(1)-1);
score_gradient = max(sum(vec_weight.*vec_filter)/(length(vec_weight)-1),0);
score = score_gradient*score_intensity;

去掉低score的Corner后，最后剩余的Corner就比较干净了，如下图。

这些角点可能还会有干扰点，这在下一步生成棋盘中会慢慢剔除。下一篇介绍另外一个重要的函数：chessboardsFromCorners。

参考资料

1、Geiger A, Moosmann F, Car Ö, et al. Automatic camera and range sensor calibration using a single shot[C]//Robotics and Automation (ICRA), 2012 IEEE International Conference on. IEEE, 2012: 3936-3943.

2、http://www.cvlibs.net/software/libcbdetect/