% RIDGEORIENT - Estimates the local orientation of ridges in a fingerprint

%

% Usage:  [orientim, reliability, coherence] = ridgeorientation(im, gradientsigma,...

%                                             blocksigma, ...

%                                             orientsmoothsigma)

%

% Arguments:  im                - A normalised input image.

%             gradientsigma     - Sigma of the derivative of Gaussian

%                                 used to compute image gradients.

%             blocksigma        - Sigma of the Gaussian weighting used to

%                                 sum the gradient moments.

%             orientsmoothsigma - Sigma of the Gaussian used to smooth

%                                 the final orientation vector field. 

%                                 Optional: if omitted it defaults to 0

% 

% Returns:    orientim          - The orientation image in radians.

%                                 Orientation values are +ve clockwise

%                                 and give the direction *along* the

%                                 ridges.

%             reliability       - Measure of the reliability of the

%                                 orientation measure.  This is a value

%                                 between 0 and 1. I think a value above

%                                 about 0.5 can be considered 'reliable'.

%                                 reliability = 1 - Imin./(Imax+.001);

%             coherence         - A measure of the degree to which the local

%                                 area is oriented.

%                                 coherence = ((Imax-Imin)./(Imax+Imin)).^2;

%

% With a fingerprint image at a 'standard' resolution of 500dpi suggested

% parameter values might be:

%

%    [orientim, reliability] = ridgeorient(im, 1, 3, 3);

%

% See also: RIDGESEGMENT, RIDGEFREQ, RIDGEFILTER

% May 2003      Original version by Raymond Thai,  

% January 2005  Reworked by Peter Kovesi           

% October 2011  Added coherence computation and orientsmoothsigma

%               made optional

% April 2018    Change computation of gradients to use DERIAVTIVE5

%

% Peter Kovesi

% peterkovesi.com

function [orientim, reliability, coherence] = ...

             ridgeorient(im, gradientsigma, blocksigma, orientsmoothsigma)

    if ~exist('orientsmoothsigma', 'var'), orientsmoothsigma = 0; end

    [rows,cols] = size(im);

    % Calculate image gradients.

    [Gx, Gy] = derivative5(gaussfilt(im, gradientsigma), 'x', 'y');

    % Estimate the local ridge orientation at each point by finding the

    % principal axis of variation in the image gradients.

    Gxx = Gx.^2;       % Covariance data for the image gradients

    Gxy = Gx.*Gy;

    Gyy = Gy.^2;

    % Now smooth the covariance data to perform a weighted summation of the

    % data.

    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);

    % Analytic solution of principal direction

    denom = sqrt(Gxy.^2 + (Gxx - Gyy).^2) + eps;

    sin2theta = Gxy./denom;            % Sine and cosine of doubled angles

    cos2theta = (Gxx-Gyy)./denom;

    if orientsmoothsigma

        sze = fix(6*orientsmoothsigma);   if ~mod(sze,2); sze = sze+1; end    

        f = fspecial('gaussian', sze, orientsmoothsigma);    

        cos2theta = filter2(f, cos2theta); % Smoothed sine and cosine of

        sin2theta = filter2(f, sin2theta); % doubled angles

    end

    orientim = pi/2 + atan2(sin2theta,cos2theta)/2;

    % Calculate 'reliability' of orientation data.  Here we calculate the area

    % moment about the orientation axis found (this will be the minimum moment)

    % and an axis perpendicular (which will be the maximum moment).  The

    % reliability measure is given by 1.0-min_moment/max_moment.  The reasoning

    % being that if the ratio of the minimum to maximum moments is close to one

    % we have little orientation information.

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

    Imax = Gyy+Gxx - Imin;

    reliability = 1 - Imin./(Imax+.001);

    coherence = ((Imax-Imin)./(Imax+Imin)).^2;

    % Finally mask reliability to exclude regions where the denominator

    % in the orientation calculation above was small.  Here I have set

    % the value to 0.001, adjust this if you feel the need

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