% Bc. Michal Kubiš
% Fakulta elektrotechniky a komunikačních technologií na VUT v Brně
% Vedoucí: Ing. Petr Fedra
% 
% 21. 5. 2010
% založeno na algoritmu Petera Kovesiho


function [orientim, reliability] = ...
             ridgeorient(im, gradientsigma, blocksigma, orientsmoothsigma)
        
    [rows,cols] = size(im);
    
    sze = fix(6*gradientsigma);   if ~mod(sze,2); sze = sze+1; end
    f = fspecial('gaussian', sze, gradientsigma); 
    [fx,fy] = gradient(f);                        
    
    Gx = filter2(fx, im); 
    Gy = filter2(fy, im); 
    
   
    Gxx = Gx.^2;       
    Gxy = Gx.*Gy;
    Gyy = Gy.^2;
    
    sze = fix(6*blocksigma);   if ~mod(sze,2); sze = sze+1; end    
    f = fspecial('gaussian', sze, blocksigma);
    Gxx = filter2(f, Gxx); 
    Gxy = 2*filter2(f, Gxy);
    Gyy = filter2(f, Gyy);
    
    denom = sqrt(Gxy.^2 + (Gxx - Gyy).^2) + eps;
    sin2theta = Gxy./denom;            
    cos2theta = (Gxx-Gyy)./denom;
       
    sze = fix(6*orientsmoothsigma);   if ~mod(sze,2); sze = sze+1; end    
    f = fspecial('gaussian', sze, orientsmoothsigma);    
    cos2theta = filter2(f, cos2theta); 
    sin2theta = filter2(f, sin2theta); 
    
    orientim = pi/2 + atan2(sin2theta,cos2theta)/2;


    Imin = (Gyy+Gxx)/2 - (Gxx-Gyy).*cos2theta/2 - Gxy.*sin2theta/2;
    Imax = Gyy+Gxx - Imin;
    
    reliability = 1 - Imin./(Imax+.001);

    reliability = reliability.*(denom>.001);