orientation.m

estimate the orientation of fingerprint

function orient= orientation(image,blocksize,noShow)

% cacluates the local flow orientation in each local window  with size (blocksize x blocksize),
% Note that the image coordinate system is x(i,vertical) axis towards bottom and y(j,horizontal) axis towards right

% input
%  image: the matrix of fingerprint
%  blocksize: the size of block 
%  noshow: don't show the figure of orientation image
% out 
% orientation: the orientation contains three element: x_label, y_label and the theta of direction

[l,w] = size(image);  % l=the number of row and   w=the number of collum for image

orient = zeros(l,w); % initialize the orientation
theta=0;

% initialize the combination of gradient 
gradient_xy = zeros(l,w);
gradient_yy_minus_xx = zeros(l,w);

% compute the y(horizontal) gradient and x(vertical) gradient through Sobel horizontal edge-emphasizing filter
filter_gradient = fspecial('sobel');
gradient_x = filter2(filter_gradient,image);   %to get x(vertical) gradient
filter_gradient = filter_gradient';
gradient_y = filter2(filter_gradient,image);     %to get y(horizontal) gradient

%compute the combination of gradient
gradient_xy = gradient_x.*gradient_y;
gradient_yy_minus_xx = (gradient_y-gradient_x).*(gradient_y+gradient_x);


xidx=1;
yidx=1;
 for i=1:blocksize:l
     xidx=xidx+1;
    for j=1:blocksize:w      
        if j+blocksize-1 < w & i+blocksize-1 < l
                            
		  sum_gradient_xy = sum(sum(gradient_xy(i:i+blocksize-1, j:j+blocksize-1)));
          sum_gradient_yy_minus_xx = sum(sum(gradient_yy_minus_xx(i:i+blocksize-1, j:j+blocksize-1)));
          
          theta=0;
          
          if  sum_gradient_yy_minus_xx ~=0 
                 theta =atan2(2*sum_gradient_xy,sum_gradient_yy_minus_xx);
                 theta=theta/2;
                 
                 orient(xidx,yidx)=theta;
                 
                 
                 
                 
                 
                 
%         
%           if  sum_gradient_yy_minus_xx>0 | sum_gradient_xy==0
%                theta = theta+pi/2;
%           end 
%            if   sum_gradient_yy_minus_xx<0  & sum_gradient_xy>0
%                  theta = theta+pi;
%            end 
          
       
        end 

        end
    sum_gradient_xy=0;
    sum_gradient_xx_minus_yy=0;  
    end 
end 

if nargin == 2
%     	imagesc(direct);
% 
% 	hold on;
	[u,v] = pol2cart(orient(:,3),blocksize);
   quiver(orient(:,2),orient(:,1),u,v,0,'g');
   hold off;
end;