关于棋盘格角点检测 Corner detection

       图像角点检测是图像视觉领域中的经典研究方向。鉴于不同类型的角点有着不同的特点,棋盘格角点作为一种特殊的角点广泛应用于摄像机标定中,角点检测在整个摄像机标定过程中也扮演着相当重要的角色。

      目前来说,主要有基于图像灰度的角点检测和图像边缘提取的角点检测,其中后一种依赖于图像边缘提取的好坏,实用性不是很强。所以当下应用最广泛的是基于图像灰度曲率变化的检测方法,常用到的算子有Moravec 、Harris,Forstner和Susan等,其中Harris的定位性能和鲁棒性较好,加之Harris算法原理简单,不受摄像机姿态和光照的影响,现已成为了使用最为广泛的角点提取算法之一。因此,针对Harris算法现在也已提出了许多改进算法:亚像素级的角点检测算法与尺度空间相结合的算法,自适应角点提取算法根据棋盘格的几何特性提出的检测算法等,但大部分都是以增加计算量或者低效率为代价,同时,用 Harris 算法进行检测,有三点不足:(1 )该算法不具有尺度不变性;(2 )该算法提取的角点是像素级的;(3 )该算法检测时间不是很令人满意。以下就介绍一种基于Harris算子的改进算法,且该算法已用到本人的论文中,效果很好。

     首先Harris 角点检测是基于图像像素灰度值变化梯度的,灰度值图像的角点附近,是其像素灰度值变化非常大的区域,其梯度也非常大。换句话说,在非角点位置邻域里,各点的像素值变化不大,甚至几乎相等,其梯度相对也比较小。如果邻域内点的灰度值与中心点Image (i,j) 的灰度值之差 的绝对值 在一个阈值t 范围 内,那就认为这个点与中心点是相似的。与此同时,属于该 Image(i,j) 点的相似点计数器nlike(i,j) 也随之加1。在 Image (i,j) 点的n 邻域 全部被遍历一边之后,就能得到在这个邻域范围内 与中心点相似的 点个数的统计值nlike(i,j) 。根据nlike(i,j) 的大小,就可以判断这个中心点是否可能为角点。由于选择3*3的检测窗口,所以,对于中心像素点 , 在下面的讨论中只考虑其8 邻域内像素点的相似度。算该范围的像素点与中心像素点的灰度值之差的绝对值 (为 Δ )  , 如果该值小于等于设定的阈值 ( 记为 t)  , 则认为该像素点与目标像素点相似。

             nlike(i,j)=sum(R(i+x,j+y))    ( -  1≤ x≤ 1 , -  1 ≤ y ≤ 1 , 且 x ≠ 0 , y ≠ 0)   ,

       其中 :  R(i+x, j+y)=1 , 当Δ( i + x , j + y)≤ t

               R(i+x, j+y)=0 , 当Δ( i + x , j + y) > t

     从定义中可以看出 : 0≤ nlike ( i , j ) ≤ 8 。 现在讨论 nlike( i , j) 值的含义 。

   (1)   nlike (i , j) = 8 ,表示当前中心像素点的 8邻域范围内都是与之相似的像素点 , 所以该像素点邻域范围内的梯度不会很大 , 因此角点检测时 , 应该排除此类像素点,不将其作为角点的候选点 。

   (2)   nlike (i , j) = 0 ,表示当前中心像素点的 8邻域范围内没有与之相似的像素点 , 所以该像素点为孤立像素点或者是噪声点 , 因此角点检测时 , 也应该排除此类像素点 。

   (3)   nlike (i , j) = 7 ,可以归结为以下的两者情况 ,其他情形都可以通过旋转来得到 ( 图中黑色区域仅表示与中心像素相似 , 而两个黑色区域像素可能是相似的 , 也可能不相似 ) 。

   (4)   nlike (i , j) = 1 ,  中心像素点也不可能为角点 。  

   (5)  2≤ nlike ( i , j) ≤ 6 , 情况比较复杂 , 无法确认像素点准确的性质。我采取的方法是先将其列入候选角点之列,对其进行计算角点响应CRF等后续操作。

   然后计算各候选角点的CBF,进行非极大抑制,再次由于Harris算子不具有尺度不变性,为确保角点检测的准确度和完整性,我们采用Mikolajczyk和Schmid提出的Harris - Laplace检测方法[8],将Harris算子与高斯尺度空间相结合,从而克服Harris 算子只能在单一尺度下检测角点的缺点

 

代码见下:

function [cim, r, c, rsubp, csubp] = myharris(im, sigma,  disp)
    
    error(nargchk(2,7,nargin));
    if nargin == 2
	   disp = 0;
    end
    
    if ~isa(im,'double')
	   im = double(im);
    end

    subpixel = nargout == 5;
    
    % Derivative masks       
    s_D = 0.7*sigma;
    x  = -round(3*s_D):round(3*s_D);
    dx = x .* exp(-x.*x/(2*s_D*s_D)) ./ (s_D*s_D*s_D*sqrt(2*pi));
    dy = dx';
    Ix = conv2(im, dx, 'same');
    Iy = conv2(im, dy, 'same');  

    % Generate Gaussian filter of size 6*sigma (+/- 3sigma) and of
    % minimum size 1x1.
    g = fspecial('gaussian',max(1,fix(6*sigma)), sigma);
    Ix2 = conv2(Ix.^2, g, 'same'); % Smoothed squared image derivatives
    Iy2 = conv2(Iy.^2, g, 'same');
    Ixy = conv2(Ix.*Iy, g, 'same');
    
    t = 25;
    Image = im;
    [nrow,ncol] = size(Image);
    Corner = zeros(nrow,ncol);
    
    %除了边界上boundary个像素,检测图像上每个点
    boundary=8;
    for i=boundary:nrow-boundary+1 
        for j=boundary:ncol-boundary+1
            nlike=0; %相似点个数
            if abs(Image(i,j)-Image(i-1,j-1))<=t 
                nlike=nlike+1;
            end
            if abs(Image(i,j)-Image(i-1,j))<=t   
                nlike=nlike+1;
            end
            if abs(Image(i,j)-Image(i-1,j+1))<=t  
                nlike=nlike+1;
            end  
            if abs(Image(i,j)-Image(i,j-1))<=t  
                nlike=nlike+1;
            end
            if abs(Image(i,j)-Image(i,j+1))<=t  
                nlike=nlike+1;
            end
            if abs(Image(i,j)-Image(i+1,j-1))<=t  
                nlike=nlike+1;
            end
            if abs(Image(i,j)-Image(i+1,j))<=t  
                nlike=nlike+1;
            end
            if abs(Image(i,j)-Image(i+1,j+1))<=t  
                nlike=nlike+1;
            end

            if nlike>=3 && nlike<=6
                Corner(i,j)=1;%如果周围有0,1,7,8个相似与中心的(i,j)
                              %那(i,j)就不是角点,所以,直接忽略
            end;
        end;
    end;


    cim = zeros(nrow,ncol);
    for i = boundary:nrow-boundary+1 
      for j = boundary:ncol-boundary+1
        if Corner(i,j)==1  %只关注候选点
            M = [Ix2(i,j) Ixy(i,j);
                 Ixy(i,j) Iy2(i,j)];
         cim(i,j) = det(M)/(trace(M) + eps);
        end
      end
    end      
    
    
    [lr,lc,lmax] = findLocalMaximum(cim,3*sigma);
    thresh = 0.15*max(lmax(:));
    radius = 2*sigma;                     
   
    if nargin > 2   % We should perform nonmaximal suppression and threshold
	  if disp  % Call nonmaxsuppts to so that image is displayed
	    if subpixel
		  [r,c,rsubp,csubp] = nonmaxsuppts(cim,  radius, thresh, im);
        else
		  [r,c] = nonmaxsuppts(cim, radius, thresh, im);	
        end
	  else     % Just do the nonmaximal suppression
	    if subpixel
		[r,c,rsubp,csubp] = nonmaxsuppts(cim, radius, thresh);
	    else
		[r,c] = nonmaxsuppts(cim, radius, thresh);		
        end
      end
    end
    

   function [row,col,max_local] = findLocalMaximum(val,radius)
   
       mask = fspecial('disk',radius)>0;
       nb = sum(mask(:));
       highest = ordfilt2(val, nb, mask);
       second_highest = ordfilt2(val, nb-1, mask);
       index = highest==val & highest~=second_highest;
       max_local  = zeros(size(val));
       max_local(index) = val(index);
      [row,col] = find(index==1);
   
   


非极大抑制函数,并且进行角点坐标亚像素求精。

function [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im)

    v = version; Octave = v(1)<'5';     % Crude Octave test    
   subPixel = nargout ==4;            % We want sub-pixel locations    
    [rows,cols] = size(cim);
    bordermask = zeros(size(cim));
    bordermask(radius+1:end-radius, radius+1:end-radius) = 1;
    
    % Find maxima, threshold, and apply bordermask
    cimmx = (cim==mx) & (cim>thresh) & bordermask
        
    for i=8:rows-7
     for j=8:cols-7
         if cimmx(i,j)==1;
             for x=i-7:i+7
                 for y=j-7:j+7
                     if cimmx(x,y)==1
                         if cim(x,y)>cim(i,j)
                             cimmx(i,j)=0;
                         end
                     end
                 end
             end
         end
     end
end
    
    
       [r,c] = find(cimmx);                % Find row,col coords. 
    
    
    if subPixel        % Compute local maxima to sub pixel accuracy  
	if ~isempty(r) % ...if we have some ponts to work with
	
	ind = sub2ind(size(cim),r,c);   % 1D indices of feature points
	w = 1;         % Width that we look out on each side of the feature
                       % point to fit a local parabola
	
	% Indices of points above, below, left and right of feature point
	indrminus1 = max(ind-w,1);
	indrplus1  = min(ind+w,rows*cols);
	indcminus1 = max(ind-w*rows,1);
	indcplus1  = min(ind+w*rows,rows*cols);
	
	% Solve for quadratic down rows
	cy = cim(ind);
	ay = (cim(indrminus1) + cim(indrplus1))/2 - cy;
	by = ay + cy - cim(indrminus1);
	rowshift = -w*by./(2*ay);       % Maxima of quadradic
	
	% Solve for quadratic across columns	
	cx = cim(ind);
	ax = (cim(indcminus1) + cim(indcplus1))/2 - cx;
	bx = ax + cx - cim(indcminus1);    
	colshift = -w*bx./(2*ax);       % Maxima of quadradic

	rsubp = r+rowshift;  % Add subpixel corrections to original row
	csubp = c+colshift;  % and column coords.
	else
        rsubp = []; csubp = [];
	end
    end
    
    if nargin==4 & ~isempty(r)     % Overlay corners on supplied image.
	if Octave
	    warning('Only able to display points under Octave');
	    if subPixel
		plot(csubp,rsubp,'r+'), title('corners detected');
	    else
		plot(c,r,'r+'), title('corners detected');
	    end
	    [rows,cols] = size(cim);	    
	    axis([1 cols 1 rows]); axis('equal'); axis('ij')
	else
	    figure, imshow(im,[]), hold on
	    if subPixel
     		plot(csubp,rsubp,'r+'), title('corners detected');
	    else	    
		plot(c,r,'r+'), title('corners detected');
	    end
	end
    end

实验结果如下:

   

         

                                                               

 

 


 

 

 

 

 

 

 

 

 

 

 

 

 

没有更多推荐了,返回首页