• 相机One Shot标定

    1 原理说明

    原理部分网上其他文章[1][2]也已经说的比较明白了，这里不再赘述。

    2 总体流程

    参考论文作者开源的Matlab代码[3]和github上的C++代码[4]进行说明（不得不说还是Matlab代码更优雅）

    论文方法总体分两部，第一部是在画面中找到所有的类棋盘格角点，第二步是角点的基础上构建出棋盘格形状。

    3 模块说明

    3.1 寻找角点

    论文寻找角点的思想是用2x4个模板对输入图像进行逐像素匹配，得到类棋盘格角点，然后使用优化方法对角点进行亚像素优化得到最终角点；

    寻找角点算法流程如下：

    总体来说寻找角点的流程可分为4步：

    1. 模板匹配：包括选择尺度、模板创建、模板匹配、结果融合
    2. 非极大值抑制
    3. 角点优化：包括角点方向计算和优化、角点位置优化
    4. 分数计算和角点调整：包括分数计算、角点调整

    3.1.1 模板匹配

    论文中模板匹配使用了2x4个模板（见上述）对全图像进行匹配，为了适应不同尺度的棋盘格图像，论文创建了3种不同尺度的2x4个模板。

    1. % 3 scales
    2. radius(1) = 4;
    3. radius(2) = 8;
    4. radius(3) = 12;
    5.  
    6. % template properties
    7. template_props = [
    8.   0 pi/2 radius(1); pi/4 -pi/4 radius(1);     % 小尺度
    9.   0 pi/2 radius(2); pi/4 -pi/4 radius(2);     % 中尺度
    10.   0 pi/2 radius(3); pi/4 -pi/4 radius(3)];    % 大尺度
    11.  
    12. disp('Filtering ...');
    13.  
    14. % filter image
    15. img_corners = zeros(size(img,1),size(img,2));
    16. for template_class=1:size(template_props,1)
    17.   
    18.   % create correlation template
    19.   template = createCorrelationPatch(template_props(template_class,1),template_props(template_class,2),template_props(template_class,3));
    20.   
    21.   % filter image according with current template
    22.   img_corners_a1 = conv2(img,template.a1,'same');
    23.   img_corners_a2 = conv2(img,template.a2,'same');
    24.   img_corners_b1 = conv2(img,template.b1,'same');
    25.   img_corners_b2 = conv2(img,template.b2,'same');
    26. end

    在得到不同模板的结果后需要进行结果融合，论文给的公式是：

    对应的实现是：

    1.   % compute mean
    2.   img_corners_mu = (img_corners_a1+img_corners_a2+img_corners_b1+img_corners_b2)/4;
    3.   
    4.   % case 1: a=white, b=black
    5.   img_corners_a = min(img_corners_a1-img_corners_mu,img_corners_a2-img_corners_mu);
    6.   img_corners_b = min(img_corners_mu-img_corners_b1,img_corners_mu-img_corners_b2);
    7.   img_corners_1 = min(img_corners_a,img_corners_b);
    8.   
    9.   % case 2: b=white, a=black
    10.   img_corners_a = min(img_corners_mu-img_corners_a1,img_corners_mu-img_corners_a2);
    11.   img_corners_b = min(img_corners_b1-img_corners_mu,img_corners_b2-img_corners_mu);
    12.   img_corners_2 = min(img_corners_a,img_corners_b);
    13.   
    14.   % update corner map
    15.   img_corners = max(img_corners,img_corners_1);
    16.   img_corners = max(img_corners,img_corners_2);

    3.1.2 非极大值抑制

    这个就不详细介绍了，反正就是非极大值抑制...

    3.1.3 角点优化

    论文中对于角点优化说的不咋详细，所以主要参考作者的开源代码。这一部分主要的步骤是：

    1. Sobel算子计算图像x,y方向梯度
    2. 根据图像梯度计算图像梯度方向和权重（强度）
    3. 根据图像附近区域像素梯度方向和权重，计算角点方向
      1. 对应论文中的梯度方向直方图统计过程
    1. function [v1,v2] = edgeOrientations(img_angle,img_weight)
    2.  
    3. % init v1 and v2
    4. v1 = [0 0];
    5. v2 = [0 0];
    6.  
    7. % number of bins (histogram parameter)
    8. bin_num = 32;
    9.  
    10. % convert images to vectors
    11. vec_angle  = img_angle(:);
    12. vec_weight = img_weight(:);
    13.  
    14. % convert angles from normals to directions
    15. vec_angle = vec_angle+pi/2;
    16. vec_angle(vec_angle>pi) = vec_angle(vec_angle>pi)-pi;
    17.  
    18. % create histogram
    19. angle_hist = zeros(1,bin_num);
    20. for i=1:length(vec_angle)
    21.   bin = max(min(floor(vec_angle(i)/(pi/bin_num)),bin_num-1),0)+1;
    22.   angle_hist(bin) = angle_hist(bin)+vec_weight(i);
    23. end
    24.  
    25. % find modes of smoothed histogram
    26. [modes,angle_hist_smoothed] = findModesMeanShift(angle_hist,1);
    27.  
    28. % if only one or no mode => return invalid corner
    29. if size(modes,1)<=1
    30.   return;
    31. end
    32.  
    33. % compute orientation at modes
    34. modes(:,3) = (modes(:,1)-1)*pi/bin_num;
    35.  
    36. % extract 2 strongest modes and sort by angle
    37. modes = modes(1:2,:);
    38. [foo idx] = sort(modes(:,3),1,'ascend');
    39. modes = modes(idx,:);
    40.  
    41. % compute angle between modes
    42. delta_angle = min(modes(2,3)-modes(1,3),modes(1,3)+pi-modes(2,3));
    43.  
    44. % if angle too small => return invalid corner
    45. if delta_angle<=0.3
    46.   return;
    47. end
    48.  
    49. % set statistics: orientations
    50. v1 = [cos(modes(1,3)) sin(modes(1,3))];
    51. v2 = [cos(modes(2,3)) sin(modes(2,3))];
    1. 角点方向优化
      1. 步骤3中计算角点方向比较粗糙，这里对方向进行优化处理
      2. 优化方法对应论文公式(4)，主要是计算角点附件有效点的梯度累加矩阵的最小特征值对应的特征向量（有点拗口:D）
    1.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2.   % corner orientation refinement %
    3.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4.   
    5.   A1 = zeros(2,2);
    6.   A2 = zeros(2,2);
    7.   
    8.   for u=max(cu-r,1):min(cu+r,width)
    9.     for v=max(cv-r,1):min(cv+r,height)
    10.       
    11.       % pixel orientation vector
    12.       o = [img_du(v,u) img_dv(v,u)];
    13.       if norm(o)<0.1
    14.         continue;
    15.       end
    16.       o = o/norm(o);
    17.       
    18.       % robust refinement of orientation 1
    19.       if abs(o*v1')<0.25 % inlier?
    20.         A1(1,:) = A1(1,:) + img_du(v,u) * [img_du(v,u) img_dv(v,u)];
    21.         A1(2,:) = A1(2,:) + img_dv(v,u) * [img_du(v,u) img_dv(v,u)];
    22.       end
    23.       
    24.       % robust refinement of orientation 2
    25.       if abs(o*v2')<0.25 % inlier?
    26.         A2(1,:) = A2(1,:) + img_du(v,u) * [img_du(v,u) img_dv(v,u)];
    27.         A2(2,:) = A2(2,:) + img_dv(v,u) * [img_du(v,u) img_dv(v,u)];
    28.       end
    29.       
    30.     end
    31.   end
    32.   
    33.   % set new corner orientation
    34.   [v1,foo1] = eig(A1); v1 = v1(:,1)'; corners.v1(i,:) = v1;
    35.   [v2,foo2] = eig(A2); v2 = v2(:,1)'; corners.v2(i,:) = v2;
    1. 角点位置优化
      1. 优化方法对应论文公式(2)，求解方法对应论文公式(3)
    1.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2.   %  corner location refinement  %
    3.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4.   
    5.   G = zeros(2,2);
    6.   b = zeros(2,1);
    7.   for u=max(cu-r,1):min(cu+r,width)
    8.     for v=max(cv-r,1):min(cv+r,height)
    9.       
    10.       % pixel orientation vector
    11.       o = [img_du(v,u) img_dv(v,u)];
    12.       if norm(o)<0.1
    13.         continue;
    14.       end
    15.       o = o/norm(o);
    16.       
    17.       % robust subpixel corner estimation
    18.       if u~=cu || v~=cv % do not consider center pixel
    19.         
    20.         % compute rel. position of pixel and distance to vectors
    21.         w  = [u v]-[cu cv];
    22.         d1 = norm(w-w*v1'*v1);
    23.         d2 = norm(w-w*v2'*v2);
    24.         
    25.         % if pixel corresponds with either of the vectors / directions
    26.         if d1<3 && abs(o*v1')<0.25 || d2<3 && abs(o*v2')<0.25
    27.           du = img_du(v,u);
    28.           dv = img_dv(v,u);
    29.           H = [du dv]'*[du dv];
    30.           
    31.           G = G + H;
    32.           b = b + H*[u v]';
    33.         end
    34.       end
    35.     end
    36.   end
    37.  
    38.   % set new corner location if G has full rank
    39.   if rank(G)==2
    40.     corner_pos_old = corners.p(i,:);
    41.     corner_pos_new = (G\b)';
    42.     corners.p(i,:) = corner_pos_new;
    43.     
    44.     % set corner to invalid, if position update is very large
    45.     if norm(corner_pos_new-corner_pos_old)>=4
    46.       corners.v1(i,:) = [0 0];
    47.       corners.v2(i,:) = [0 0];
    48.     end
    49.     
    50.   % otherwise: set corner to invalid
    51.   else
    52.     corners.v1(i,:) = [0 0];
    53.     corners.v2(i,:) = [0 0];
    54.   end

    3.1.4 分数计算和角点调整

    • 分数计算

    分数计算是用来计算每个角点的得分，然后与用户设置的阈值进行比较，如果超过阈值则为合格角点。这里论文说的和代码里面的有点对不上感觉，论文说方法的是求角点附近图像的二阶导，然后和一个模板做normalized cross-correlation；但是代码里面实现的方式是：

    1. 计算论文中Figure 2.e中的模板T
    1. % center
    2. c = ones(1,2)*(size(img_weight,1)+1)/2;
    3.  
    4. % compute gradient filter kernel (bandwith = 3 px)
    5. img_filter = -1*ones(size(img_weight,1),size(img_weight,2));
    6. for x=1:size(img_weight,2)
    7.   for y=1:size(img_weight,1)
    8.     p1 = [x y]-c;
    9.     p2 = p1*v1'*v1;
    10.     p3 = p1*v2'*v2;
    11.     if norm(p1-p2)<=1.5 || norm(p1-p3)<=1.5
    12.       img_filter(y,x) = +1;
    13.     end
    14.   end
    15. end
    1. 将模板T与之前得到的图像梯度权重分别做归一化，然后做相关，得到相关系数score_gradient
    1. % convert into vectors
    2. vec_weight = img_weight(:);
    3. vec_filter = img_filter(:);
    4.  
    5. % normalize
    6. vec_weight = (vec_weight-mean(vec_weight))/std(vec_weight);
    7. vec_filter = (vec_filter-mean(vec_filter))/std(vec_filter);
    8.  
    9. % compute gradient score
    10. score_gradient = max(sum(vec_weight.*vec_filter)/(length(vec_weight)-1),0);
    1. 按角点的2个主方向生成4个模板
    2. 计算模板匹配值score_intensity（过程与3.1.1的模板匹配一样）
    1. % create intensity filter kernel
    2. template = createCorrelationPatch(atan2(v1(2),v1(1)),atan2(v2(2),v2(1)),c(1)-1);
    3.  
    4. % checkerboard responses
    5. a1 = sum(template.a1(:).*img(:));
    6. a2 = sum(template.a2(:).*img(:));
    7. b1 = sum(template.b1(:).*img(:));
    8. b2 = sum(template.b2(:).*img(:));
    9.  
    10. % mean
    11. mu = (a1+a2+b1+b2)/4;
    12.  
    13. % case 1: a=white, b=black
    14. score_a = min(a1-mu,a2-mu);
    15. score_b = min(mu-b1,mu-b2);
    16. score_1 = min(score_a,score_b);
    17.  
    18. % case 2: b=white, a=black
    19. score_a = min(mu-a1,mu-a2);
    20. score_b = min(b1-mu,b2-mu);
    21. score_2 = min(score_a,score_b);
    22.  
    23. % intensity score: max. of the 2 cases
    24. score_intensity = max(max(score_1,score_2),0);
    1. 计算score = score_intensity * score_gradient
    2. 改变模板尺度（与 3.1.1一样），重复步骤1~5，计算所有尺度中的分数最大值

    不知道论文和代码中的方法是否等价，留以后考证。

    • 角度调整

    这部分主要是对计算角点的主方向做一个调整，以减少后面多图像匹配（我们用不到）过程中的歧义性

    1. % make v1(:,1)+v1(:,2) positive (=> comparable to c++ code)
    2. idx = corners.v1(:,1)+corners.v1(:,2)<0;
    3. corners.v1(idx,:) = -corners.v1(idx,:);
    4.  
    5. % make all coordinate systems right-handed (reduces matching ambiguities from 8 to 4)
    6. corners_n1 = [corners.v1(:,2) -corners.v1(:,1)];
    7. flip       = -sign(corners_n1(:,1).*corners.v2(:,1)+corners_n1(:,2).*corners.v2(:,2));
    8. corners.v2 = corners.v2.*(flip*ones(1,2));

    3.2 构建棋盘格

    论文提出的棋盘格构建是基于生长的方法，所以首先要选定几个点作为初始棋盘格，然后在初始棋盘格的基础上向外生长，在生长的过程中需要对多种可能的生长方式进行判断，判断新生长出来的棋盘格是不是最优，最终得到最优的棋盘格结构。

    通过上述步骤，可以得到多个棋盘格结构，如果这些棋盘格结构中有重叠的话，则需要做一个去重处理。

    对于判断生长和去重过程中的哪个棋盘格结构更优，论文给出了一种棋盘格能量计算方法，参考论文公式(6)。

    所以整体棋盘格构建流程如下：

    3.2.1 选择初始点

    没什么讲究，所有点都拿来试一试...

    3.2.2 构建初始棋盘格

    构建初始棋盘格是将选择的初始点扩展成3x3共9个点，组成一个小小的初始棋盘格。方法是从初始点开始，往2个主方向和2个主方向的反方向（一共4个方向）各找一个距离最近的点，这样就有了5个点，然后再从上下两个点以同样的方法往左右扩展，最终形成3x3棋盘格。

    1. % return if not enough corners
    2. if size(corners.p,1)<9
    3.   chessboard = [];
    4.   return;
    5. end
    6.  
    7. % init chessboard hypothesis
    8. chessboard = zeros(3,3);
    9.  
    10. % extract feature index and orientation (central element)
    11. v1 = corners.v1(idx,:);
    12. v2 = corners.v2(idx,:);
    13. chessboard(2,2) = idx;
    14.  
    15. % find left/right/top/bottom neighbors
    16. [chessboard(2,3),dist1(1)] = directionalNeighbor(idx,+v1,chessboard,corners);
    17. [chessboard(2,1),dist1(2)] = directionalNeighbor(idx,-v1,chessboard,corners);
    18. [chessboard(3,2),dist2(1)] = directionalNeighbor(idx,+v2,chessboard,corners);
    19. [chessboard(1,2),dist2(2)] = directionalNeighbor(idx,-v2,chessboard,corners);
    20.  
    21. % find top-left/top-right/bottom-left/bottom-right neighbors
    22. [chessboard(1,1),dist2(3)] = directionalNeighbor(chessboard(2,1),-v2,chessboard,corners);
    23. [chessboard(3,1),dist2(4)] = directionalNeighbor(chessboard(2,1),+v2,chessboard,corners);
    24. [chessboard(1,3),dist2(5)] = directionalNeighbor(chessboard(2,3),-v2,chessboard,corners);
    25. [chessboard(3,3),dist2(6)] = directionalNeighbor(chessboard(2,3),+v2,chessboard,corners);
    26.  
    27. % initialization must be homogenously distributed
    28. if any(isinf(dist1)) || any(isinf(dist2)) || ...
    29.    std(dist1)/mean(dist1)>0.3 || std(dist2)/mean(dist2)>0.3
    30.   chessboard = [];
    31.   return;
    32. end

    可以看到，在通过directionalNeighbor扩展了之后，还对扩展出去的距离做了一个判断，避免产生太过离谱的初始棋盘格。

    directionalNeighbor的具体实现如下，其中dist_point+5*dist_edge作为距离这一个有点迷，可能是想兼顾考虑径向和切向距离，但是不知道数学原理是什么（虽然效果确实不错）

    1. function [neighbor_idx,min_dist] = directionalNeighbor(idx,v,chessboard,corners)
    2.  
    3. % list of neighboring elements, which are currently not in use
    4. unused       = 1:size(corners.p,1);
    5. used         = chessboard(chessboard~=0);
    6. unused(used) = [];
    7.  
    8. % direction and distance to unused corners
    9. dir  = corners.p(unused,:) - ones(length(unused),1)*corners.p(idx,:);
    10. dist = (dir(:,1)*v(1)+dir(:,2)*v(2));
    11.  
    12. % distances
    13. dist_edge = dir-dist*v;
    14. dist_edge = sqrt(dist_edge(:,1).^2+dist_edge(:,2).^2);
    15. dist_point = dist;
    16. dist_point(dist_point<0) = inf;
    17.  
    18. % find best neighbor
    19. [min_dist,min_idx] = min(dist_point+5*dist_edge);
    20. neighbor_idx = unused(min_idx);

    3.2.3 棋盘格生长

    这个流程没啥好说的，确定了初始棋盘格之后，就不断往4个方向生长，然后选一个能量最低的作为最优解，作为新的初始棋盘格。

    1.   % try growing chessboard
    2.   while 1
    3.     
    4.     % compute current energy
    5.     energy = chessboardEnergy(chessboard,corners);
    6.     
    7.     % compute proposals and energies
    8.     for j=1:4
    9.       proposal{j} = growChessboard(chessboard,corners,j);
    10.       p_energy(j) = chessboardEnergy(proposal{j},corners);
    11.     end
    12.     
    13.     % find best proposal
    14.     [min_val,min_idx] = min(p_energy);
    15.     
    16.     % accept best proposal, if energy is reduced
    17.     if p_energy(min_idx)<energy
    18.       chessboard = proposal{min_idx};
    19.       
    20.     % otherwise exit loop
    21.     else
    22.       break;
    23.     end
    24.   end

    3.2.4 棋盘格去重

    当棋盘格中包含的任意一个角点在已存在的棋盘格中也存在了，则认为存在棋盘格重复。去重采用的方法是看看当前棋盘格的和已存在的棋盘格哪个更优(哪个能量更低)，如果当前棋盘格能量更低，则替换掉原来的棋盘格。

    1.   % if chessboard has low energy (corresponding to high quality)
    2.   if chessboardEnergy(chessboard,corners)<-10
    3.  
    4.     % check if new chessboard proposal overlaps with existing chessboards
    5.     overlap = zeros(length(chessboards),2);
    6.     for j=1:length(chessboards)
    7.       for k=1:length(chessboards{j}(:))
    8.         if any(chessboards{j}(k)==chessboard(:))
    9.           overlap(j,1) = 1;
    10.           overlap(j,2) = chessboardEnergy(chessboards{j},corners);
    11.           break;
    12.         end
    13.       end
    14.     end
    15.  
    16.     % add chessboard (and replace overlapping if neccessary)
    17.     if ~any(overlap(:,1))
    18.       chessboards{end+1} = chessboard;
    19.     else
    20.       idx = find(overlap(:,1)==1);
    21.       if ~any(overlap(idx,2)<=chessboardEnergy(chessboard,corners))
    22.         chessboards(idx) = [];
    23.         chessboards{end+1} = chessboard;
    24.       end
    25.     end
    26.   end

    4 参考资料

    [1] Geiger, A., Moosmann, F., Car, Ö. and Schuster, B., 2012, May. Automatic camera and range sensor calibration using a single shot. In 2012 IEEE international conference on robotics and automation (pp. 3936-3943). IEEE.
    [2] 基于生长的棋盘格角点检测方法--（1）原理介绍_findchessboardcornerssb原理介绍_计算机视觉life的博客-CSDN博客
    [3] Andreas Geiger (cvlibs.net)
    [4] onlyliucat/Multi-chessboard-Corner-extraction-detection-: chess board corner extraction and chess board recovery "Automatic Camera and Range Sensor Calibration using a single Shot" (github.com)

  • 相关阅读:
    CentOS二进制安装Containerd
    Java面向对象进阶3——多态的概述及特点
    PHP redis sorted set 有序集合
    zookeeper客户端命令
    如何让Nginx更安全？
    五、Express
    浅谈建筑能耗智能监测平台发展现状及未来趋势
    nginx.4——正向代理和反向代理（七层代理和四层代理）
    辅助解决小白遇到的电脑各种问题
    mysql json数据类型 相关函数
  • 原文地址：https://blog.csdn.net/sharemon/article/details/133148060